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"PLANS — A Finite Element Program for Nonlinear 

Analysis of Structures, Volume I - Theoretical Manual" 

by 

A. Pifko, H. S. Levine, and H. Armen, Jr. 

Grumman Aerospace Corporation 
Bethpage, New York 11714 

1. INTRODUCTION AND SUMMARY 

This report describes the current state of the PLANS system, 
a finite element program for the nonlinear analysis of structures. 
PLANS, rather than being a single comprehensive computer program, 
represents a collection of special purpose computer programs or 
modules, each associated with a distinct physical problem class. 
Utilizing this concept, each module is an independent finite ele- 
ment computer program with its associated element library that 
can be individually loaded and used to solve the problem class of 
interest. In this manner, any simplification or specialization 
germane to the individual analysis can be incorporated within each 
program thereby optimizing computational efficiency. Thus, un- 
avoidable inefficiencies that result from a more general purpose 
program are avoided in each special purpose module. 

All the programs in PLANS employ the "initial strain" concept 
within an incremental procedure to account for the effect of plas- 
ticity and include the capability for cyclic plastic analysis. 
Geometric nonlinearities included in several of the modules are 
treated by using the "updated" or convected coordinate approach. 
These nonlinear terms may be incorporated as effective loads 
and/or modifications to the stiffness matrix. 



The increased demands from governmental and industrial or- 
ganizations associated with aerospace, naval, and nuclear reactor 
technology for determining accurate stress, strain, and displace- 
ment fields have been a motivating force behind the activity and 
advances in the development of techniques for the nonlinear anal- 
ysis of structures. 

As a result of these advances, the computational capability 
available for the nonlinear analysis of structures has experi- 
enced a tremendous growth during the past ten years. Indeed, the 
level of structural analysis capability that was achieved has out- 
stripped our ability to describe accurately complex material be- 
havior such as cyclic-, time-, and temperature -dependent plas- 
ticity. Prior to the development of the programs now available, 
the designer or analyst, confronted with a problem involving 
structural nonlinearities, was left with a choice of using his 
engineering judgment alone, or in conjunction with potentially 
expensive laboratory tests. He now has the further option of 
performing numerical analysis to gain insight into the behavior 
of the structure. 

Most nonlinear analysis programs, with the exception of a 
few, have been developed as a spin-off of existing programs that 
were originally designed for linear structural analysis. Although 
this development is a natural one, the added dimension of general- 
ity has placed great responsibilities on the user of such pro- 
grams. Perhaps the greatest asset of these programs, the ability 
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to solve sophisticated problems, also represents a potential lia- 
bility, i.e., they always produce numbers. The user must still 
exercise engineering judgment in order to interpret the results 
meaningfully. The analytic results will confirm these feelings 
and provide him with additional insight. However, he now has a 
luxury of having his intuition fail him without suffering the 
consequences of a catastrophic failure, or an overdesigned system. 

The progress associated with solving nonlinear problems has 
been the direct result of advances in several interdisciplinary 
areas that include structural mechanics, numerical analysis, and 
computer systems engineering. Specifically, in the area of struc- 
tural mechanics, significant advances have been those associated 
with the finite element method. The period of these advances for 
linear elastic behavior can be traced to the development of the 
displacement method of finite element analysis (Ref. 1). 

References 2-16 are representative of advances associated 
with incorporating the effects of material nonlinear behavior 
into a finite element analysis. The approach taken in these 
references is to develop solutions of the repetitive type that 
linearize the basically nonlinear problem for a sequence of 
loading steps. These developed techniques fall into two cate- 
gories: the initial strain (Refs. 2-10) and tangent modulus 

(Refs. 11-16) methods. They differ computationally depending 
on whether the effect of material nonlinearities enters as 
pseudoloads (initial strain) or the stiffness matrix is ex- 
plicitly altered (tangent modulus) . 

Considerable effort has also been directed towards the treat- 
ment of geometric nonlinearities. These efforts are reported in 
Refs. 17-22 for consideration of geometric nonlinearities alone 
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while the simultaneous treatment of both types of nonlinearities 
was reported in Refs. 23 - 21 . 

The concurrent contribution from numerical analysis methods 
notably has been the development of efficient algorithms to solve, 
on a repetitive basis, large systems of equations. The ease with 
which the theoretical considerations and numerical algorithms are 
efficiently processed and translated into meaningful answers is 
attributable to the ingenuity and talents of computer systems en- 
gineers who develop the necessary hardware and, in some cases, 
software to satisfy what must seem to them the insatiable appe- 
tite of the developers of large scale programs. 

The advances in the above-mentioned areas have resulted in 
the development of many software systems to treat nonlinear be- 
havior. These range from simple instructional programs to gen- 
eral comprehensive computer systems capable of treating the en- 
tire visible spectrum of two and three dimensional problems. 
Reference 28 presents in some detail the current capability of 
the various general purpose finite element programs that are 
available for plastic analysis while Ref. 29 details the avail- 
able geometric nonlinear programs. 

This report describes the theoretical basis of the PLANS sys- 
tem, a collection of special-purpose finite element computer pro- 
grams for nonlinear analysis. These programs are an outgrowth of 
the work cited in Refs. 6, 8, and 30 conducted by the authors on 
the development and implementation of finite element methods for 
nonlinear analysis. As such, in order to give an historical pro- 
spective of the program's evolution, we first chronologically 
summarize the development of the capabilities represented by 
PLANS. 
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1.1 Chronology of the Development of Nonlinear Structural 
A nalysis Techniques at Grumman Aerospace Corporation 


Contract; 

Period; 

Ref. Document; 
Goal; 

Result: 

Comments : 


NAS 1-5040 
5/65-12/66 
Ref. 6 

To combine existing finite element tech- 
nology with plasticity theories to treat 
general (including cyclic) loading con- 
ditions . 

Methodology developed for membrane 
stressed structures (two dimensional). 
Some verification with NASA experiments 
on notched panels under cyclic loading. 

The plastic analysis program developed 
at this stage required a separate elastic 
analysis program to generate influence 
coefficients. A restriction of a maximum 
of 55 elements with approximately 100 
degrees of freedom was enforced on the 
basis of core storage limitations. This 
early effort is notable for its intro- 
duction of a reasonably sophisticated 
plasticity theory (kinematic hardening) 
within the framework of a complex numeri- 
cal structural analysis procedure. The 



Contract : 
Period ; 

Ref. Document: 
Goals: 


Results : 


Comments : 


Contract : 
Period: 


sjmthesis of this capability has served 
as a forerunner for subsequent efforts. 

NAS 1-7315 
6/67-6/69 
Ref. 8 

1) Extend methods developed in previous 
study to treat nonlinear bending be- 
havior; and 2) determine the feasibility 
of treating combined material and geo- 
metric nonlinearities. 

1) Introduced and applied the concept of 
an elastic-plastic boundary in the plane 
and through the thickness of triangular 
and rectangular bending elements. 

Treated problems associated with bending 
alone or combined with membrane loads. 

2) Incorporated, on a limited basis, an 
incremental technique to account for 
combined material and geometric non- 
linearities . 

Applications to beams and plates. De- 
termined collapse loads and modes of a 
variety of plate-type structures. Com- 
bined nonlinearities created for beams 
and arches. 

NAS 1-10087 
6/70-6/74 



Ref. Document: 


Ref. 30 and current report 

Goals: 1) Expand "library” of finite elements to 

include a broad base of applications; 

2) develop general comprehensive program 
for the plastic analysis of structures; 

3) extend the methods for the treatment 
of combined material and geometric non- 
linearity; and 4) implement combined 
nonlinear analysis techniques into the 
program developed in (2). 

Results: 1) In addition to triangular membrane and 

plate elements incorporated in the first 
tvjo studies, the following elements were 
included into the element library: axi- 

symmetric revolved triangle, axisymmetric 
shell, 3-D variable node isoparametric 
solid element, shear panel, stringers, 
axisymmetric rings, general beams of vari- 
ous cross sections, and a composite ele- 
ment. 2) PLANS has evolved as a collec- 
tion of programs each one of which is de- 
signed to treat a particular class of 
analysis, i.e., axisymmetric analysis of 
bodies of revolution, 3-D analysis of 
solid bodies, etc. (the particular mod- 
ules are discussed in more detail in a 
subsequent section of this report). 

3) General techniques to treat combined 
nonlinearities have been developed within 
the framework of the previous procedures 
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to treat plasticity alone, 4) The com- 
bined nonlinear analysis techniques were 
implemented into PLANS for the analysis 
of axisymmetric shell structures and for 
general 3-D built-up structures, typical 
of aircraft construction. 

Comments: The successful completion of the goals 

associated with this effort has enabled 
analysts to realistically consider prob- 
lem areas that could otherwise be only 
treated on a heuristic basis. Several 
of these problem areas are discussed sub- 
sequently in this section. 

1 . 2 Applicable Problem Areas — Current and Future 

Only two related motivating factors could justify the cost 
associated with developing advanced analytic capabilities such as 
those represented by PLANS. The first motivation is that the in- 
vestment will result in reducing the cost associated with design- 
ing and/or building the system. The second reason for developing 
advanced analyses is that the system being considered does not 
fall into a general category of those previously designed and in 
service. Therefore, there must be the capability to gather in- 
sight by test and analysis to demonstrate the integrity of the 
system. 

In addition to the problems associated with reducing weight 
and treating overload conditions, the following problem areas ap- 
pear to be the most likely to benefit — from the viewpoint of the 
above-mentioned motivating factors — by the further development 
of nonlinear structural analyses. 
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1.2.1 Fatigue and Fracture Analysis 

Fatigue life predictions are by-and-large semi empirical pro- 
cedures. The techniques generally used have, until recently, 
paid little more than lip service to the effects of plastic flow. 
It is not realistic to consider such phenomena as the effect of 
residual stresses, propagating cracks, and cyclic load spectrum 
in a flawed structure without a reasonably accurate treatment of 
plastic deformation. Also included in the fatigue and fracture 
analysis problem area is fastener technology, where the installa- 
tion and operation of the fasteners involve large strains and 
deformations . 

1.2.2 Crashworthiness Evaluation 

A capability for the analysis of general transportation sys- 
tems is valuable for evaluating existing designs, post-mortem 
analysis of damaged units, predicting the crashworthiness of pro- 
posed designs, and establishing crashworthiness design criteria. 
Results of crash simulation studies can be used to determine prob- 
able damage to passengers, equipment, and structure. In addition, 
areas for structural modifications (including energy absorbers) 
can be determined along with their associated cost and/or weight 
penalties. Note that the application of such a capability is not 
limited to transportation systems, but include design and analysis 
of such stationary structures as highway barriers and nuclear con- 
tainment vessels subjected to impact loads. 

1.2.3 High Temperature Applications 

In the design and analysis of nuclear reactor system compo- 
nents, large thermal gradients are common; in fact, much of the 
plastic deformation results from thermally induced loadings. 
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Furthermore, these loadings tend to be cyclical, separated by 
hold periods under sustained high levels of constant stress or 
strain conditions. The space shuttle represents another example 
of a structure that vjill operate in an environment in which tem- 
peratures are sufficiently high to cause degradation of struc- 
tural strength. For both the reactor components and the space 
shuttle it becomes difficult to separate the effects associated 
with plastic flow and creep. 

1.2.4 Metal Forming Technology 

Ductile metals may be rolled, hammered, bent, or extruded 
even at relatively low temperatures. These processes generally 
require considerable energy. Therefore accurate methods for de- 
termining the behavior of solids, under the conditions severe 
enough to cause permanent changes in shape, are significant in 
attaining increased fabrication efficiency. 

1.2.5 Nonisrotropic Materials 

The treatment of the nonlinear behavior of structures con- 
structed of homogeneous orthotropic or layered composite materi- 
als may appear to be presumptuous in view of the many remaining 
unanswered questions associated with homogeneous isotropic mate- 
rials. Nevertheless such materials exist (more so than the ideal- 
ized materials considered) and are being used with increasing fre- 
quency in aerospace applications. It would be unrealistic to ne- 
glect the special problems posed by such materials when developing 
a general purpose plastic analysis capability. 
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2. PROGRAM DEFINITION 


2 . 1 General Discussion 

Developing a comprehensive finite element program remains a 
subjective undertaking since it invariably depends on the analysts 
designing the program, the computer hardware available, and the 
resources allocated for the project. 

Many of the basic criteria for both linear and nonlinear 
general purpose codes are similar. However, given that the ap- 
propriate theories from structural mechanics have been imple- 
mented, the distinguishing features of a nonlinear program are 
twofold: 1) displacements, stresses, elastic-plastic strains, 

etc., must be stored for subsequent use in succeeding load steps, 
and 2) the solution is of the repetitive type so that calcula- 
tions that are performed once for an elastic analysis must be re- 
peated in a nonlinear analysis. The latter consideration is most 
critical because while it is true that nonlinear analyses have 
reached the stage of maturity such that almost any problem can be 
solved for a price, the cost for such an analysis can become pro- 
hibitive. Consequently, it is Incumbent on the designer of a 
nonlinear code to minimize the cost of an analysis so that solu- 
tions to meaningful problems are economically feasible. 

For a large comprehensive program, steps towards this end 
must be taken not only in reducing computing time (i.e., use of 
the central processing unit) but also in reducing overhead costs 
associated with accessing data on secondary storage devices. 

The use of an efficient solution algorithm has the most sig- 
nificant effect in reducing computer costs associated with CPU 
time . 
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After the choice of a solution algorithm has been made, great 
care must be taken in coding the key "number crunching" portion of 
the code, including the judicious use of machine language routines 
for certain operations. Efforts must also be made to perform key 
input/output functions in an optimum fashion so that their asso- 
ciated overhead costs are minimized. The allocation bet\«?een 
scratch disc storage and primary in-core storage is the main op- 
erational decision to be made here. Such features as variable 
core allocation enable this ratio to be optimized for a given 
problem and machine. 

2 . 2 The PLANS System 

The PLANS system, rather than being one comprehensive com- 
puter program, is actually a collection of finite element programs 
used for the nonlinear analysis of structures. This collection of 
programs evolved and is based on the organizational philosophy in 
which classes of analysis are treated individually based on the 
physical problem class to be analyzed. On the basis of this con- 
cept, each of the independent finite element computer programs of 
PLANS, with an associated element library can be individually 
loaded and used to solve the problem class of interest. A number 
of programs have been developed for material nonlinear behavior 
alone and for combined geometric and material nonlinear behavior. 
Table 1 summarizes the usage, capabilities, and element libraries 
of the current programs of the PLANS system. These include: 

• REVBY for the plastic analysis of bodies 

of revolution 

• OUT-OF-PLANE for the plastic analysis of built-up 

structures where membrane effects 
predominate 
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Table 1 


CURRTNT CAPABILITY OF PLANS 


Module 

Usage 

Element Library 

Loading 

Current Capablllcv 

OUT -OF -PLANE 
(material nonllnearicy) 

Analysis of built-up 
structures composed of 
membrane skins, 
sf-ingers, bulkheads 

Family of triangular membrane 
elements, stringer element, 
shear panel, beam elements, 
with the following cross sec- 
tions: solid rectangular, 

solid circular, I-section, 
Z-section, L-section, 
T-section, hollow rectangular, 
hollow circular 

- Concentrated forces and 
moments at nodes 

- Distributed edge load on 
triangular element 

- Distributed lateral load on 
beam elements 

- Thermal loading 

600 members 
900 nodes 

6 degrees of freedom 
per node 

kEVBY 

Miaterial nunlineari ry) 

Anaiyois of bodies of 
revolution 

Isoparametric shell of revolu- 
tion element, revolved tri- 
angular element, ring stiffener 
of various cross sections 

Ax i symmetric 

- Line loading 

- Distributed load on shell 
elements 

- Distributed load on revolved 
triangles 

- Thermal loading 

600 members 
900 nodes 

6 degrees of freedom 
per node 

SATELLITE 

Mesh generation, plot 
j Healizarion , data 
checking 


- 

■ 

0- PLANE -MG 

(maceLial and geometric 
nonllneafi ty) 

Same as OUT -OF -PLANE but 
includes the effect of 
geometric nonlinear 
behavior 

Same as OUT -OF -PLANE 

Same as OUT -OF -PLANE 

Same as OUT-OF- PLANE 

HEX 

(material nonlinearity) 

Analysis of three dimen- 
sional solids 

Isoparametric family of 
hexahedra 

- Concentrated forces at nodes 

- Distributed surface load on 
a face of a hexahedra 

- Thermal loading 

2300 members 
2500 nodes 

3 degrees of freedom 
per node 

BEND 

(material I’Onl .iviearity) 

Analysis of built-up 
structures, includes 
bending of sheet material 

Family of triangular membrane 
elements, stringer element, 
beam element (same cross sec- 
tions as OUT -OF- PLANE ) , higher 
order triangular plate element 
(bending and membrane) 

- Concentrated forces at nodes 

- Distributed edge loads on 
triangular place and mem- 
brane elements 

- Distributed lateral load on 
plate and beam elements 

- Thermal load 

Open core - variable num* 
ber of members and nodes 
Up to 12 degrees of free- 
dom per node 



BEND 


for the plastic analysis of built-up 
structures where bending and membrane 
effects are significant 

• HEX for problems requiring a three dimen- 

sional elastic, plastic analysis 

• OUT -OF -PLANE -MG for the material and geometric non- 

linear analysis of built-up structures 

Supplementing these is a SATELLITE program for data debugging and 
plotting of input geometries. 

Developing the programs in this manner afforded some distinct 
advantages, particularly in a research environment, since we were 
able to keep the individual programs as simple as possible and, 
therefore, easily understood and modifiable. Furthermore, in 
general, each program follows the same basic flow and each con- 
tains many common subroutines. Consequently, once the basic flow 
and supervisory subroutines have been established, development 
could proceed as parallel efforts with different programs being 
worked on at the same time . 

It also simultaneously allowed for the development of spin- 
off modules for special purpose analysis. An example of this, re- 
ported in Ref. 31, is the FAST module for fracture analysis. An- 
other special purpose module developed under this concept is 
SATELLITE which integrates the input part of the programs with 
FORTRAN mesh generation subprograms and input plotting routines. 

Another implication of this approach also arose in the de- 
velopment of PLANS. With the development of each succeeding pro- 
gram, new and hopefully better programming techniques were real- 
ized. These improvements were implemented in succeeding programs 
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so that continued growth in program capability, generality, and 
efficiency was achieved. 

Although the programs can be individually loaded it is con- 
venient to integrate them into a single system by making each 
module callable from an executive program. 

A brief discussion of the significant features of several 
components of the system follows. 

2.2.1 Executive Program 

The executive program is a small calling program that uti- 
lizes the computer's system to form an overlay structure. How- 
ever, for an IBM facility we make use of a PL/1 program called 
WIZARD (Ref. 32), written at Grumman Data Systems, that allows 
a user to store many program decks on a magnetic tape or direct 
access volume. 

This system is comprised of three separate collections of 
subprograms. These are: 

• A source update program to allow a user to selectively 
edit and compile a new source program 

• An object update program that maintains updated files 
of compiled source programs 

• A file select program that selectively loads programs 
from the object file for execution 

While the file select program can accept commands to load 
individual programs, it can also load groups of programs by simply 
supplying a group name that has been previously defined. In addi- 
tion, each group definition may contain other group names in its 
definition. This feature is particularly meaningful for our 
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purposes since each program naturally defines a group, and pro- 
grams common to all groups (such as the solution package) are 
naturally subgroups. 

The modules of the PLANS system are also operational on CDC 
computers. However, since PL/1 is unavailable on many CDC facili- 
ties, the system described above cannot be used. Consequently, on 
CDC machines we have made use of an overlay approach to implement 
the executive system. Specifically, use is made of the NASTRAN 
linkage editor facility to implement this overlay procedure. 

2.2.2 Basic Flow of the Analysis Programs 

Modules were developed in the PLANS system that treat mate- 
rial nonlinear behavior and combined material and geometric non- 
linear behavior. 

The modules that analyze material nonlinearities alone im- 
plement the initial strain approach which does not require the 
stiffness matrix be updated at any step in the analysis so that 
the program organization will differ from that for combined geo- 
metric and material nonlinear behavior. First, we discuss in 
detail the basic flow for a typical module in PLANS that treats 
material nonlinear behavior alone using the initial strain ap- 
proach. 

Figure 1 is a schematic representation of the flow of each of 
the analysis programs that implement this approach. It should be 
pointed out that the program outlined below does not fulfill the 
criterion of being general purpose in the sense that it can be 
used to perform analysis of different physical phenomena (geo- 
metric and material nonlinear behavior) and implement various solu- 
tion procedures (tangent modulus and initial strain approach). 
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Program in PLANS 
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However, it was developed with the intent to implement one solu- 
tion procedure for a given problem class in the simplest and most 
cost effective manner. In this sense the entire PLANS system and 
not an individual module can be considered to be a general pur- 
pose program. 

As shown in Fig. 1, each module has three major components: 
a main calling routine (MAIN), an elastic analysis subprogram 
(ELAS) , and the plastic analysis subprogram (PLAS) . A discussion 
of each of these components follows. 

• MAIN 

Each individual computer program is controlled by a main 
calling program that sets up core allocation for principal arrays 
and the data set specifications for auxiliary storage. The main 
program then transfers control to subroutine ELAS. ELAS is a 
finite element elastic analysis program that performs the elastic 
analysis and calculates the initial yield load. Control is trans 
ferred back to the main program, which calls subroutine PLAS only 
if requested. PLAS manages the plastic analysis and maintains 
control of the analysis until the complete plastic analysis is 
performed. By so organizing the program it should be possible 
with a minimum amount of change to use PLAS with other available 
finite element elastic analysis programs. 

e ELAS 

Figure 1 shows a block diagram of the computational flow of 
ELAS. This program is a special purpose finite element program 
for the elastic analysis of structures. Accordingly, its first 
major task is to read all input. The input is read in functional 
groups as follows: 
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• Problem title 

• Nodal coordinates and control variables 

• Member topology describing connectivity 

• Nodal boundary conditions. Single and multipoint con- 

straints on the displacements are specified in addi- 
tion to fixed or free conditions . 

• Load vector — includes member or nodal load data for 
distributed surface loads, edge or line loads, concen- 
trated loads , and thermal loads 

• Material and section properties — tables of elastic and 
plastic material properties and member geometric proper- 
ties are set up along with applicable members. Complete 
generality has been maintained so that orthotropic mate- 
rial behavior can be specified. 

Other variables controlling output are also specified. The input 
scheme has been written so as to minimize the amount of card 
handling for any problem. 

During the input portion of the program, externally numbered 
nodes and members are converted to internally numbered nodes and 
members via some internally generated tables. As a consequence 
of this feature, nodes and members can be arbitrarily numbered and 
new ones can be easily added to an existing idealization or the 
connectivity (bandwidth) altered with the specification of some 
simple card input. In addition, based on these tables, there is 
a substantial amount of input checking for the validity of speci- 
fied nodes and members. 

The next step is to form all element stiffness, stress, and 
initial strain stiffness matrices. This routine also places the 
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elements of the stiffness matrix along with their position in the 
total stiffness matrix on a sequential auxiliary storage device. 
In addition, the element stress and initial strain matrix along 
with some element control variables are placed on another sequen- 
tial auxiliary storage device to be used subsequently for stress- 
strain calculations. 

The total stiffness matrix and load vector are assembled by 
sequentially reading the element stiffness matrix and "stacking" 
indices from auxiliary storage and stacking that portion of the 
stiffness matrix that fits in core, and then reading it onto an 
auxiliary storage device. This process is repeated untij. the en- 
tire load vector and stiffness matrix have been formed. 

Our basic design philosophy is to make alterations, perform 
certain matrix manipulations such as those required for single 
and multipoint constraints (see Sec. 3.6), and account for bound- 
ary conditions (including applied displacements) on the element 
level, and then to assemble these quantities into the master ar- 
rays. In this manner we need not make use of a matrix interpre- 
tive system to manipulate large matrices. 

With the total stiffness matrix and load vector assembled, 
the next step is to solve for displacements. Two solution pack- 
ages are currently being used for this purpose, both based on the 
Cholesky factorization method. Their differences are based on 
the system used to store the stiffness matrix. PIRATE is based 
on partitioning the structure so that the stiffness matrix is ex- 
plicitly in block tridiagonal form. The second routine, PODSYM, 
is based explicitly on the banded nature of the stiffness matrix, 
so that only elements within the lower triangle that lie between 
the semibandwidth need be assembled. In addition, this routine 
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"packs” the rows of the stiffness matrix before writing it on 
auxiliary storage by suppressing all consecutive zeros. These 
algorithms are presented in detail in Sec. 3.7 of this manual. 

The unit load stresses are calculated next from displacements. 
These stresses are used to calculate the initial yield load, and 
this yield load is used to scale the displacements and calculate 
yield load stresses and strains. Control is then returned to the 
main program. 

• PLAS 

Figure 2 shows the computational flow of a typical subroutine 
PLAS. This program supervises the entire plastic analysis after 
initial processing has been carried out by ELAS. The principal 
information that is communicated to PLAS is: 

• Factored stiffness matrix and unit load vector 

® Element stress and initial strain matrices 

® Initial yield load 

• Plastic material properties 

With this information, the load is incremented and the re- 
petitive calculations that implement the incremental solution are 
performed as follows: 

• Add incremental effective plastic load vector including 
any contribution from equilibrium correction factor to 
scaled vector of applied loads 

• Solve for increments in displacements using the fac- 
tored stiffness matrix. This involves only a forward 
and backward solution of a triangular matrix. 
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Fig. 2 Flow of 


Subroutine SICREM supervises 
reading sequential auxiliary 
storage devices containinf, 
element stresses, total 
strains, etc. and calculates 
Increments in stress, strain, 
and other pertinent quanti- 
ties. The plasticity consti- 
tutive relations are imple- 
mented. The effective plastic 
load vector is formed, includ- 
ing the equilibrium correc- 
tion factor. Total quantities 
are formed by summing incre- 
mental quantities. Print 
stresses, strains if re- 
quested . 


Calculate "zero" load 
stresses and strains and dis- 
placements. Find new critical 
load and calculate new critical 
load stresses and strains. 


i.ne PLAS 








• Read total displacements from disc storage and add dis- 
placement increments to form updated total displace- 
ments. Write results on disc storage. Read total ele- 
ment (or nodal) stress, strains, plastic strains (shift 
in the yield surface when kinematic hardening is used) 
from disc storage. Calculate element (or nodal) incre- 
ments of total strain, stress, etc., check yield cri- 
terion for previously elastic elements (nodes) and un- 
loading condition for elements currently in the plastic 
range, implement plastic constitutive relations and form 
effective plastic load vector. Sum total quantities 
with calculated incremental quantities and write on disc 
storage to be used as input for the next Increment of 
load . 

® Repeat the above if the maximum load has not been 
reached. 

When the maximum load is reached, the final effective plastic 
load vector is formed and saved. At this point, a check is made 
to see if another cycle of loading is desired. If not, control 
is returned to the main program. If an additional cycle of load- 
ing is specified, new plastic material parameters can be read as 
input and inserted into the material property tables. Displace- 
ments are then calculated with the applied load vector set equal 
to zero. Stresses and strains are calculated on the basis of 
these displacements. If subsequent yielding occurs in the re- 
versed cycle at a load that is opposite in sign to the load pre- 
viously computed, then these stresses and strains represent re- 
sidual quantities. This is checked next by computing a new yield 
load. This calculation involves the solution of a quadratic 
equation for each member. One root of this solution represents 
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the previously reached maximum load if no new plastic material 
properties were input, and the other represents the new yield 
load in the reversed direction. This load level is used as a 
starting point for the next cycle of loading. The effective 
plastic load vector is then added to the applied load vector, 
and the sum is used to solve for the new critical load displace- 
ments and the corresponding stresses and strains. Transfer is 
then made to the beginning of PLAS in order to increment the load 
and proceed as described previously. 

2 , 3 Combined Material and Geometric Nonlinear Analysis 

For combined material and geometric nonlinear analysis sever- 
al major alterations in the program flow must be made because the 
convected coordinate approach to the solution of geometric non- 
linear problems requires the reformation of the stiffness matrix 
during the incremental solution process. This flow is presented 
in Fig. 3. As shown in this illustration, no distinction can be 
made between an ELAS and a PLAS as before. The main program con- 
trols the flow previously supervised by ELAS and PLAS. Now there 
is a special routine to initialize stresses, strains, etc., where- 
as previously all quantities were initialized to be the yield load 
stresses, strains, and displacements. Aside from this distinction 
the flow of both analyses is similar. Now for the geometric non- 
linear analysis, after an increment of load has been applied, in- 
crements of displacement are calculated and the geometry is up- 
dated. Subroutine SRAIN is then called. This subroutine super- 
vises the elastic-plastic calculations. In addition to calculating 
the element stresses, strains, etc., using the appropriate stress- 
strain relations, the element stiffness matrices and mechanical 
load vector are updated because of the -geometry changes and the 
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Fig. 3 Flow of Program for Combined Material and 
Geometric Nonlinear Behavior 
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presence of initial stresses. These quantities are stored on 
auxiliary storage devices and control is returned to the main 
program. The new incremental load vector, including plastic and 
geometric effects as well as the "equilibrium correction," is 
then fo'rmed. If the maximum load has not been reached, control 
returns to point A or B in the flow depending upon whether 
the stiffness matrix is to be updated and refactored. If the 
maximum load for this load sequence has been reached new plastic 
material properties (and possibly a new load vector) are read and 
the new load increment is calculated. For reversed loading there 
can be no elastic scaling to the next yield load since the geome- 
try changes make the unloading a nonlinear problem. Consequently 
the load sense is reversed, and the load is incremented to the new 
minimum load. 


3. THEORETICAL APPROACH 

This section describes the theoretical basis of the programs 
of the PLANS system. Included in this section are the formula- 
tion of the governing matrix equations that account for material 
and geometric nonlinear behavior, the description of the plas- 
ticity theory implemented in the programs, and some computational 
aspects of the programs. Much of the material in this section is 
explained in great detail in the authors' previous NASA contractor 
reports (Refs. 6, 8, and 30), and so the discussion in this sec- 
tion has been restricted to a review of the highlights. The 
above-mentioned reports may be referred to for additional details. 
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3 . 1 Governing Matrix Equations 


The approach used by the authors in several investigations, 
spanning a number of years, for the plastic analysis of structures 
(Refs. 6, 8, 9, 10, 26, and 30), incorporates the initial strain 
concept for the treatment of material nonlinear behavior and the 
incremental moving coordinate formulation for the treatment of 
geometric nonlinear behavior. As such, this is the approach im- 
plemented in the programs of the PLANS system. An extensive pre- < 
sentation of the governing equations for this approach is pre- 
sented in Ref. 8 and for completeness is reviewed below. 

If the increments in displacements (AU], linear total 

strains (Ae..}, and rotations (Aoj..} are related to the nodal 
i-J ij 

generalized displacements following matrix rela- 

tions 

(AU) = [N](AU^} 

{Ae_} = [W]{AU^} (1) 

{Acd,.} = [a]{AU^] 

and is some initial stress state, then the matrix equation 

governing the nonlinear response of a structural finite element 
to some arbitrary increment of loading can be written as (Ref. 26) 

([k®] + [k^l)CAU^} = CAP^} + CAQ^) + (R^) (2) 

where 

[k^] = conventional stiffness matrix obtained from the 

linear component of the strain-displacement rela- 
tions 
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geometric stiffness matrix, formulated from the 
nonlinear terms of the strain-displacement rela- 
tions. It can be considered as an additional 
component of the stiffness matrix that accounts 
for the effect that the presence of initial 
stresses have on subsequent deformations. 

increments of generalized nodal forces 

increments of the effective plastic load. If an 
assumption is made on the plastic strain distribu- 
tion v?ithin the element then 

initial strain stiffness matrix used to account 
for the presence of initial strains, and reflects 
the assumed distribution of both the total and 
plastic strains within the element 

increments of nodal or element plastic strain 

vector of residual forces that may exist because 
of equilibrium imbalance at the end of the pre- 
vious load increment. 


[W]^[E] [W]dV 



V 


(AQo) 


tl 


[W]‘^[E]{Ae^}dV 


Equation (2) is derived by using a moving coordinate system 
fixed to the body. It is valid for large elastic-plastic deforma- 
tions, provided the appropriate nonlinear terms are retained in 
the strain-displacement relations and the total strain increment 
can be simply decomposed into elastic and plastic components. 
Additionally, proper transformations from the previous to current 
coordinate systems must be used so that the changes in orienta- 
tion and volume of the elements are accounted for (Refs. 26 and 
24). In its usage herein, consideration is restricted to small 
strain, moderate rotation problems. 

3 .2 Solution Procedures 

3.2.1 Small Deflection Problems (Material Nonlinearity Alone) 

The matrix formulation presented in Eq. (2) is used for this 
case with the [k ] matrix deleted. The individual element 
stiffness matrices and load vectors are assembled to form the 
global stiffness matrix. Boundary conditions and single and 
multipoint constraint conditions are recognized at this stage and 
accounted for at the element level, thereby removing any further 
consideration of these conditions on the larger assembled coef- 
ficient matrices. At this stage we have 

= (AP^}^'''^ + {AQ^}’- + {R^}^ (3) 

Here the superscripts i and i+1 refer to the current and 
next load step. Since the exact values of the plastic strain in- 
crements for the next load step are not known, we make use of 
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values from the current step. The equilibrium correction is 
based on total equilibrium at the end of the current load step 
and is used in every increment. Note that the stiffness matrix 
need never be reformed after the first increment since the non- 
linearities appear as a component to the load vector. Conse- 
quently, the stiffness matrix can be factored once for the ini- 
tial load step and used repeatedly throughout the incremental 
solution. A significant reduction in solution time (factors of 
3 to 21) have been realized depending on the size of the 
problem and the bandwidth using the previously factored coeffi- 
cient matrix. The solution procedure is as follows. 

A specified "unit" load is applied to the structure and cor- 
responding stresses and strains are determined. Since the re- 
sponse is linear up to initial yield, a critical load for which 
plastic deformation first occurs at a node (or element) can be 
calculated from the unit load stresses. From this level the load 
is incremented to a maximum value with new increments of dis- 
placement, plastic strain and stress calculated at each step. 
Total values are obtained by summing incremental values. 

If the load is reversed for a cyclic analysis, new material 
properties data may be read in at this point. A new critical 
load for which yielding begins in the reverse direction is cal- 
culated, based on elastic unloading to this point. Procedures 
for determining this load and initial yield load are presented 
in Refs. 6 and 8. The critical load for reversed loading may 
occur before all the initial load is removed because of the pre- 
sence of residual stresses and the existence of the Bauschinger 
effect. The load is then incremented from the new critical value 
to a specified maximum (minimum) value. This procedure is re- 
peated for as many half cycles as desired. 
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3.2.2 Combined Material and Geometric Nonlinearity 

Again the formulation presented in Eq. (2) is employed. The 
load is applied in small increments from the initial unloaded 
state. At the end of each increment, new increments of deflec- 
tion, stress, total strain, and plastic strain are calculated. 
Total quantities are summed and appropriate transformations from 
previous to current geometry are used to calculate such quanti- 
ties as the initial stresses, geometry of the struc- 

ture is updated in each step. Again the plastic strain incre- 
ments used for the next step are those calculated at the end of 
the current step. 

Now the element contributions are assembled and the system 
of linear incremental equations is solved for the next load incre- 
ment. For the large deflection problem, the most time-consuming 
feature is the reassembly of the stiffness matrix in every incre- 
ment, and the necessity of resolving the equations. To minimize 
this time-consuming feature it becomes convenient to treat the 
large deflection terms as effective loads and not reform the 
stiffness matrix [k^] in every increment (Ref. 26). We can then 
rewrite Eq. (2) as 

= - [k^HAU^}’- + CAP^)^'’’^ + (AQ^}^ + {R^]^ (4) 

Now the product of the geometric stiffness times the current dis- 
placement increment is treated as an "effective load." This is 
less accurate than the "tangent modulus" formulation but the 
"equilibrium correction" prevents excessive drifting of the solu- 
tion. The stiffness matrix [k^] is updated every M steps 
(M > 1) where M is input by the user. Additionally, if large 
non linearities or instability are anticipated by the user, he may 
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switch to the tangent modulus approach for geometric nonlinearity 
at any point in the analysis. After this specified load, the 
stiffness matrix [k] = [k^] + [k^] is reformed and resolved for 
every increment of load. 

The solution process is repeated until the maximum specified 
load is reached or structural failure occurs. If cyclic loading 
is specified, the load increment is reversed at the maximum load 
(a new magnitude may be input) , new material properties data are 
read, as well as new values of M and the crossover load for 
tangent modulus treatment of geometric nonlinearities. The in- 
cremental process is then repeated until the new maximum (mini- 
mum) load is reached and the procedure is repeated for as many 
load cycles as desired. 

3 .3 Plasticity Relations 

This section considers appropriate incremental plasticity re- 
lations to determine values of stress and plastic strain developed 
during the history of loading. In PLANS use is made of Hill's 
yield criterion (Ref. 33), for an orthotropic material, which re- 
duces to the von Mises yield condition for isotropic materials, to 
predict initial yield and obtain the flow rules of plasticity. 

The capability of handling both strain hardening and ideally 
plastic behavior is included in the program. 

3.3.1 Matrix Relations - Strain Hardening 

We can, for small strain increments, decompose the total 

T e 

strain increment {Ae }' into elastic {Ae } and plastic {Ae} 
components, as 

(Ae"^} = (Ae®} + (Ae] (5) 
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The elastic strain increments are related to the stress incre- 
ments (Aa) by 

{Ae®} = [E]”^[Aa} (6) 

where [E] ^ is an array whose elements are combinations of 
elastic material constants. 

A linear incremental constitutive relation between plastic 
strains and stresses can be written from all of the popularly 
used flow theories of plasticity. This relationship can be 
represented as 

[As] = [C](Aa3 (7) 

Therefore, substituting Eqs. (6) and (7) into Eq. (5) we can write 

(Aa} = [R]'^(Ae'^} (8) 

where 

[R] = [E]’^ + [C] 

The elements of the matrix [C] depend on the current state 
of stress, yield condition, flow rule, and hardening law. Exam- 
ples of the explicit form of the [C] matrix are given in 
Tables 2 and 3. In PLANS the yield condition is based on Hill's 
yield criterion for an orthotropic material (Ref. 33), and the 
hardening law is based on the Prager-Ziegler kinematic hardening 
theory (Refs. 34, 35, and 36) modified for orthotropic material 
behavior. Both linear and nonlinear strain hardening options 
are available with input parameters determining which is chosen. 

To minimize input requirements for nonlinear hardening, a Ramberg- 
Osgood (Ref. 37) representation of the stress-strain data is 
used 
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(9) 


e 


a , 3a f a 

E 7E Va 


n-1 


0.7 


Thus, for this representation of the stress-strain la\j, for 
an initially isotropic material two additional material parame- 
ters n and Oq y are required for the plastic analysis. For a 
Ramberg-Osgood representation of an initially orthotropic mate- 
rial these two quantities are required for each principal mate- 
rial direction and for each of the shear components. For a three 
dimensional problem this means that six sets (three for the nor- 
mal components and three for the shear components) of stress- 
strain data must be specified. The flow rule is based on 
Drucker's postulate (Ref. 38). 


3.3.2 Matrix Relations - Perfect Plasticity 

The treatment of multiaxial elastic ideally plastic behavior 
required that the following conditions be satisfied: 

® The stress increment vector must be tangent to the 
loading surface 


• The plastic strain increment vector must be normal to 
the loading surface, where the loading surface is the 
representation in stress space of the initial yield 
function or the subsequent yield function after some 
plastic deformation has occurred. 

If f(a^j) represents the yield surface, the first condition 
can be expressed analytically as 




( 10 ) 
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Table 2 Isotropic iCl Matrix for Various Stress States 


PLANE STRESS 


[C] 


symmetric 


m, — p * 2 ^ 


m« = a 

2 y 


= 3t 

3 xy 


n 3 2 

D = 2 ==^0 


yield function: 


Oq « yield stress , 

P - 2^-2 - - ^ ,-2 2 . 

z - o +a -oa +3 t -o =0 
X y X y xy o 


unloading criterion: m^da + m„da + m^dx < 0 

^ 1 X 2 y 3 X- 


xy 


THREE DIMENSIONAL 


[C) = 






m.m„ rn. m 


symmetric 


‘4 3 


m^m, m^mn m-m^ m-m, 

51 52 53 54 


62 63 64 65 


1 X y z 

2 y ■ ^ 2 x^ 

m <3 = o - 7 l(r + n ) 

3 z ‘ X y' 

m. = 3 t 

4 xy 

= 3t 

5 xz 


m^ = 3x 
6 yz 


yield function: 


D = :r CO , o = yield stress 
2 o o 


,2 


f = 


, ,-2 , ,-2 , ,-2 2 - 

^ — t: + 75 + s ^ + 3 t + 3t + 3 r - m = 0 

2 2 2 xy xz yz o 


unloading criterion: m^do + nudo + m< 3 do + m, dx + ni_dx + m.dx < 0 

^ Ixzy3z4xy5xz6yz 


ONE NORMAL STRESS - TWO SHEAR COMPONENTS 
2 


[Cl 


symmetric 

2 


m« = 3x 

2 xy 

m^ = 3x 

3 xz 


D = 2 > *^0 ~ yield stress , o^. = 


ij 

-2 -2 -2 2 
f = o + 3x^ + 3x -0=0 

X xy xz o 


yield function: 

unloading criteria: m,do + m«dx + m«dx < 0 

lx 2 ^y 3 xz 
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Table 3 Orthotropic 


[C] Matrix for Various Stress States 


PLANE STRESS 


[C] = i 




symmetric 

2 


m^m^ m2ni2 


2 

2 . ™3 


m. = 2(G+H)a - 2Ha 

1 X y 

m. = 2(F+H)a - 2Ho 

2 y X 


m^ = 4 Nt 

3 xy 

m, = -2Fa - 2Go 

4 y X 


D = c(mj^ + m2 + — + m^) , 

—9 

yield function; 
unloading criterion: "*■ *^2^^y ^S'^^xy ^ 


'ij 


f = (G4-H)oJ + (F+H)0^ - 2Ho^0y + 2NxJy = 1 


THREE DIMENSIONAL 


[C] 


™1™2 


m^m^ ^2^3 


ni 4 m , m , 

14 2 4 


m^m^ 


LL-, Ul^ lU^UL^ 

16 2 6 


rn~ m . 
3 4 


m2tn^ m^m^ 


3 6 


symmetric 


m, m_ 
4 5 


m-m^ m«m^ m„m^ 


4 6 5 6 


m / iti — 

A , 5 . 6 n 


2 2 ^ t+ ^ vj 

D = c(m^ + m2+m^ + -^ + -2^ + *^) , 


m. = 2(G+H)o - 2Ha - 2Ga 

1 ^ ^ X y z 

m. = 2(F+H)a - 2Fo - 2Ha 

2 ^ ^ y z X 

m^ = 2(F-FG)a - 2Go - 2Fa 

3 ^ ^ z X y 


m, = 4 Lt 
4 yz 


m c = 4 Mt 
5 zx 


m, = 4 Nt 
6 xy 


'ij °ij ■ "ij 

yield function: 

unloading criterion: + ^2^^y ^3^*^z "^4‘^^yz ^6^^xy ^ 

CH-H = 1/X^ , H+F = 1/Y^ , F-rC = 1/Z^ 

2L = 1/R^ , 2M - 1/S^ , 2N = 1/T^ 

X>Y,Z are yield stresses in tension in principal directions 

R,S,T are yield stresses in shear in principal directions 


f = F(0y -0^)2 + G(0^ - 0^)2 + H(0^ - 0y)2 + 2LXy^ + 2Mx^^ + 2 n7^^ = 1 
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and provides a linear relationship among the components of stress 
increment. Thus, one of the components tnay be expressed in terms 
of the others. In matrix form, this can be written as 

CAcr) » [EHM} (11) 

where represents the independent stress components. 

The normality condition provides a linear relation among the 
various components of the plastic strain increment. This condi- 
tion is derived from the flow rule and provides a linear rela- 
tionship in which each of the components of plastic-strain incre- 
ment can be written in terms of any one component. This rela- 
tionship may be represented in the following form 

{Ae} = [E]{A€] (12) 

f\j r\u 

where {Ac} is the independent plastic strain increment. 

The independent increments of stress and plastic strain can 
be combined and written as the components of a vector, {Aou} 

(see Ref. 8), so that Eqs . (11) and (12) can be written, respec- 
tively, as 


{Aa} = [EHA03} 
{Ae} = [E]{A(X>} 


(13) 


Examples of the explicit form of [E] and [E] are given in 

I'w 

Tables 4 and 5. Combining the above equations with Eqs. (5) and 
(6) , we can form the following relation for the independent 
quantities 

-1 

{Aco} = [E*] (Ae'^} (14) 


where 

[E*] = [E]'^[E] + [E] 
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Table 4 Isotropic [E] and fE] for Various Stress States 
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. [E] = 

0 

0 



X 
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THREE DIMENSIONAL 



"’l “ ^ ’ '"2 " '^z ^1 

m3 = +ay))/,-n^ , = ST^^/tn^ 


= 3t / m ^ 
5 xz I 


, m. = 3 t /m., 
6 yz 1 


ONE NORMAL STRESS “ TWO SHEAR COMPONENTS 



0 

-mi 




1 

0 

O' 

= 

0 

1 

0 

, [E] 

= 

mi 

0 

0 


0 

0 

1 



'"2 

0 

0 


mi 

= 3t 

xy 

/^x 

f 

= 3t 

yz 
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PLANE STRESS 



0 

-mi 


1 

0 

0 

[E] = 

0 

1 

0 , [E] = 

mi 

0 

0 


0 

0 

1 

^2 

0 

0 


m, = ( (G+H)a - Ha ' /( (F+H)o - Ha 

1 y x/ \ X y 

m„ = 2L-r / ; (F+H)o - Ha : 

2 xy \ X y/ 


THREE DIMENSIONAL 




G+H = 1/x^ 

H+F = 1/y^ 

, F+G = 1/Z^ 


2L = I/R^ 

2M = 1/S^ 

, 2N = 1/T^ 

X,Y,Z 

are yield stresses 

in tension 

in principal directions 

R,S,T 

are yield stresses 

in shear in 

principal directions 
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These relations are used with Hillfe yield criterion for ini- 
tially orthotropic behavior which reduces to the von Mises yield 
criterion for initially isotropic materials . 

3 . 4 Hardening Coefficient 

One of the unresolved questions in the use of the kinematic 
hardening theory of plasticity (Refs. 34, 35, and 36) is the 
definition of the hardening coefficient, c, (appearing in [C], 
Tables 2 and 3) under multiaxial stress states. This is espe- 
cially true for problems involving nonlinear hardening under 
cyclic loading. Several definitions have been proposed by the 
current authors (Refs. 6 and 8) and successfully used for ini- 
tially isotropic material behavior. An extension of this defi- 
nition for initially orthotropic behavior based on Hill's yield 
criterion has been recently proposed (Ref. 39) . 

Before discussing the hardening coefficient, let us first 
rewrite the basic equations required for a plastic analysis using 
kinematic hardening. 

First, we need a yield condition 

= 0 (15) 

Hill's yield criterion for orthotropic materials (Ref. 33), 
which reduces to the von Mises yield condition for isotropic 
materials, is used in PLANS 

f = F(a - a )^ + G(a^ - o )^ + H(o - a )^ (16) 

\ y g/ ' z X ' X y ^ 

+ 2Lt^ + 2M?^ + 2NT^ = 1 

yz zx xy 

where 


40 



G+H 




~ = H+F 2M = — 

Y S 

= F+G 2N = ^ 

Z T 


Here a., = a.. - a., with a., being the components of the 
shift in the yield surface and X,Y,Z are the yield stresses in 
tension in the principal directions of orthotropy and R,S,T are 
the yield stresses in shear with respect to the principal axes. 
For initially isotropic materials this reduces to the familiar 
von Mises yield condition 


f = 


- 2 

a - a 




- 2 
a “ cf . 

V z ^, 




a - cf - 
z x^ 


2 2—2 2 
+ 3t + 3t + 3t - <y =0 
xy yz zx o 


(17) 


with 0 ^ equal to the initial yield stress in tension. 

Stability requirements for use of Hill's form of the yield 
surface, i.e., ensuring that it remains a "closed ellipsoidal" 
surface, require that the following criteria must also be satis- 
fied by the yield stresses (Ref. 40) 


F F 

^ir 22 


Fi2 > 0 


^ 22^33 


^23 > ° 


F F 
^ 33^11 


> 0 


(18) 


where 
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F 


11 
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If Eq. (18) is not satisfied by the yield stresses then some other 
yield criterion must be chosen (see Ref. 40, for example) and imple 
mented in the code. This occurs frequently for composite materials 

The flow rule used in PLANS is based on Drucker's postulate 
(Ref. 38) 



dA 


do. . 


(19) 


where dA is a positive scalar quantity. For kinematic hardening 
it is given by 



df 


da 


da 


mn 


mn 


df 




kl 


■)( 


Vdaki^' 


( 20 ) 


and c is the hardening coefficient to be defined. The incre- 
mental shift in the yield surface assumed to be along 
a radial line through the center of the yield surface as postu- 
lated by Ziegler (Ref. 36) 


da . . = du (a. . - a . . ) 
ij ^ ij iJ 


( 21 ) 


and 
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dy, = 


( 22 ) 


Bf 


Ba. 




da. . 
3-J 


/ ^ 
fa. . - a . . ; r — 

^ ±2 iJ Ba 


ij 


In evaluating the terms of Eqs . (19) through (22) for sub- 
sets of the nine dimensional stress space, all terms should be 
retained in the yield function. Symmetry and setting stresses 
identically to zero should be done after the expressions have 
been obtained. The resulting expressions are termed complete 
kinematic hardening. 


3.4.1 Kinematic Hardening Coefficient for Initially Isotropic 
Materials 


If v?e evaluate c for a uniaxial state of stress and re- 
strict ourselves to initially isotropic materials, using Eqs. (19), 
(20), and (17), we get 


P 

c 2 da 


(23) 


If we generalize to a multiaxial state of stress, we may 

P ^ ~P 

consider a and e to be effective quantities a and e 

[Eq. (9)] where for an isotropic material 


~2 

a 




, cr 

/ _x 


, cr - a. 

, 

-r \ o / 


a " a , 
z x^ 


9 2 ? 

+ 3t + 3t + 3t 
yz zx xy 


(24) 


and 




P P 
2/3 de . .d€ . . 
ij iJ 
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mil I 



a 


Note the difference in definition of the "effective stresses," 
and the yield function of Eq. (17). 


If we use a Raniberg -Osgood representation of the stress- 
strain curve and consider the nonlinear term to represent the 
plastic strain contribution, then 


1 

c 


3 de^ 


da 


3 

2 




[n-l^ 

3n 

Co 

a 


7E 

Q 

O 

• 

i 


or for linear strain hardening 


(25) 


1 

c 


3 de 


da 


3 

2 


1 

IE 


1 - 


E/ 


Ej/E 


(26) 


where E^ is the tangent modulus. 

Various methods of choosing an appropriate value for E^/E 
from a stress-strain curve have been suggested (Ref. 41, for ex- 
ample) . The use of linear strain hardening can lead to substan- 
tial underpredication of plastic strains, in the range of strain 
where the material does not have a linear stress-plastic strain 
relationship. In turn, this can lead to serious discrepancies 
if the response to cyclic loading is desired. 


To determine the values of n and Oq ^ that best fit the 
actual stress-strain data, the method suggested by Ramberg- Osgood 
can be used if the strain range is sufficiently small. For this 
case 


n , 1 + l°s(17/7) 

The quantities Oq -j and Oq are 

the curve has secant moduli of 0.7E and 


) 


the stresses at which 
0.85E, respectively. 
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If the strain of interest is sufficiently large so that the pa- 
rameters as determined by the preceding process do not fit the 
curve well, then a power law representation to fit the actual 
data can be used 


€ = I + Pa” (27) 

Now having determined p and n to "best fit" the experi- 
mental data, the value of n input is that calculated, and the 
value of Oq y input is derived by equating the nonlinear terms 
of Eqs. (27) and (9), i.e.. 


1 



For cyclic loading the problem of how to determine c for 
load reversals subsequent to the initial loading is even more 
vexing especially for nonlinear hardening. For linear hardening 
a procedure for selecting the value of c and new yield stresses 
for succeeding cycles is suggested in Ref. 41. For nonlinear 
hardening it is desirable to reproduce the actual hysteretic 
stress-strain curves obtained from test data. Morrow (Ref. 42) 
details extensive cyclic tests on isotropic metals to determine 
the cyclic and hysteresis stress-strain curves. He has found 
that both the cyclic and hysteresis curves may be well repre- 
sented by a power law similar to the Ramberg-Osgood curves. 

Based on his findings, we have incorporated into PLANS a harden- 
ing coefficient for load reversals following the initial mono- 
tonic loading defined as 
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1 3 d€" 3 3ii a 


^ 2^^ 2 7E 


where we now define a to be 


* -k 
o - a V 


k k 

a - a . 


a_ - or 


2J + V 2 J 


z\ . / Z X 


V 2 


2 2 2 

„ 

+ 3t + 3t + 3t 
xy yz zx 


vjhere a^j = , and 0 ^^ is the last value of stress at 

the end of the previous loading range. 


3.4.2 


Orthotropic Hardening Coefficient 


A generalization of the definition of c for initially 
orthotropic materials based on Hill's yield criterion is pre- 
sented in Ref. 39. The basic highlights viill be reproduced here 

We define c to be 


Sf _ 

, T C>0 , 

1^ _ X 

c c 2f c 
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xy yz zx 


|(G+H)2 + h2 + g 2| d€^ 


(29) 


(F+H)2 ^ 

(F+H)2 + h2 + f 2| d€y 
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c = 
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' ^ 
9 9 9^ P 

(F+G) + F + G d€ 

j z 


c =2 

xy 




c - 2 


dT 

— O - 

. P 
ay 
^xy 

5 

yz 

9 

zx , 

dy 


This form reduces to the individual uniaxial stress-strain la^ws 
and the initially isotropic case given in Eq. (23). 

For a linear strain hardening approximation we use 
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where Em , E,p , etc., are specified by the user for each com- 

X ^y 

ponent direction. For a Ramberg -Osgood form of the uniaxial 
stress-strain data we use 
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de' 

5 
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3n 
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xy 
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(31) 

(Cont . ) 




(F+G+H) ~2f 
2L 


0.7 


yz 




0.7 


zx 


n -1 
zx 


vihere 


-2 
cr = 


2 (F+G+H) 


F(a - 0 )^ + G(0 - 0 )^ + H(0 ~ 0 ) 

^ y z' ^ z X X 


(32) 

2 2 2 1 

+ 2Lt + 2Mt + 2Nt 

yz zx xy J 

is the effective stress for initially orthotropic materials. Now 
n , 0 rt 7 5 n , 0« , etc., are specified by the user for each 

component direction. 


-2 
0 = 


For cyclic loading, we define o to be 

W l" ^^y ■ ^z^ ^ ‘^^^z ■ ^x>^ 

.2 „2 


2 (F+G+H) 


2 2 2 
^ ^ * "k -k 
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^ y z ^ z X ^ X y^ 
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vhere 


* 




A 



(34) 


Data for orthotropic cyclic loading is seriously lacking. If the 
representation given by Eqs. (33) and (34) for the effective stress 
is deemed valid then should be set to twice the measured 

0. 7ij 

monotonic value, for load reversals after the initial loading range. 
Additionally, at the end of each half-cycle of loading new plastic 
material properties may be input so that each hysteresis loop may 
be accurately represented. 


3 . 5 Unloading Criterion 

At times, due to nonproportional loading or other reasons, 
local unloading occurs even though the applied load is increasing. 
The unloading criterion, which is checked in the program for every 
load increment , is given by 


df 


da. . > 0 


UU . . / 

ij - 


^f 


oa. . xj 


da . , < 0 


for loading or neutral loading 


for unloading 


The values of used in this calculation are obtained from 

the "elastic" stress-strain relations for that increment. They 
are the actual da^j if unloading is detected. If the unloading 
criterion is not met, however, the values are determined 

from the plasticity constitutive relations. When unloading is 
detected at a point, all further stress and strain increments 
are elastic until reloading is detected using an appropriate 
yield criterion. 
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Specific relations for the unloading condition are shown in 
Tables 2 and 4. 


3.6 Single and Multipoint Constraints 


It is the philosophy of all the programs of the PLANS system 
to perform all matrix operations on the element level rather than 
on the assembled global matrices. Thus, element stress and 
strain recovery, load vector calculation, and boundary and multi- 
point constraint operations are carried out with respect to small 
matrices or vectors that readily fit in primary storage. In 
this manner, with the exception of the solution of the global 
stiffness equations, matrix operations can be accomplished with- 
out the aid of a matrix package for large scale matrix opera- 
tions . 


Accordingly, single and multipoint constraints are satisfied 
on the element level before assembling the total stiffness matrix. 
This can be illustrated by first writing the element influence 
coefficient equation as 

(£) - [k]{A 3 - [k*](€3 (35) 

o 

where 


{f] is the vector of nodal generalized forces 

{A } is the vector of nodal generalized global dis- 
® placements 

(e) is some initial strain state 

and [k], [k ] are the element stiffness and initial strain 
matrices . 
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This equation can be partitioned into an analysis set that 
contributes to the total stiffness natrix and a reaction set that 
contributes only to the calculation of nodal reaction forces. 


In PLANS, the reaction set for single point constraints 
arises from fixed node points or nodes that are given a pre- 
scribed nonzero generalized displacement component. 


Accordingly, Eq. (35) becomes 


where 



I 

k 

aa , 


ar 


k 1 k 
! rr 


ra 




Cf } = 

^ a 


[k ](A®} 
aa^'- g^ 


(36) 


(37) 


contribute to the total matrix equations with the load vector 
consisting of 

and the reaction forces are 

(f ] = [k ](A^] + [k ]{A^] - [k*](e] 
r ra g rr g l ^-ji. j 

The above procedure is accomplished within the program without 
explicitly forming the partitioned matrix by making use of an 
element location vector that is formed when the boundary condi- 
tions are initially specified. This vector is of the form 
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^ = < 
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n. 
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n . 
J 


0 


(38) 


-m, 

X 


where the zeros appear for fixed degrees of freedom, n^, n^, 
are associated v?ith degrees of freedom that are free, and -m^ 
indicates an applied displacement component. 

The use of the element location vector is demonstrated sym- 
bolically by placing the vector as shown below 
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Then all intersections of the 
dicate matrix elements that are in 


0 , 
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rr 


-m^ rows and columns in- 
and the intersection of 


0, -m. with the n's indicates matrix elements that are in k 

1 ra 


or k 

ar 

The indices 


The intersection of the n's denotes elements in k 


aa 


denote the actual row-column location in the 
total stiffness matrix of the matrix element to be assembled. In 


n . , n . 

i’ J 


52 



practice however, the total stiffness matrix is assembled in a 
one dimensional array so that n^, n^ are used to calculate a 
stacking index to relocate the matrix element in the one dimen- 
sional storage array. 

A 

The vector f is also assembled into the total load vector 
a 

using =£. In this case all elements adjacent to 0, -m. are in 

A ^ 

f and those adjacent to n, are in f and are, therefore, 

^ th ^ ^ 

added to the *'n^" location of the load vector. 

The programs of the PLANS system implement a multipoint con- 
straint capability of the form 



d t h 

where 6^ is the i dependent degree of freedom, are 

prescribed coefficients, and 6- are a set of independent 
degrees of freedom. 

Equation (39) can be used to form a transformation matrix 
which contains the pertinent a., and maps the element displace- 
ment vector {A ) which contains dependencies to a set contain- 
S 

ing only independent degrees of freedom. 

Thus , 


(A } = [T]{A 1 

O O 


(40) 


from which Eq. (35) becomes 


[f] = [k][A^3 - [k ][€} 


(41) 
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where 

Oj- m 

(f) = [T]^[f) 

[k] = [T]'^[k][T] 

[^] = [T]^[k*] 

Again the procedure is accomplished without explicitly form- 

r\j 

ing the transformation matrices or (f), [k], [k ]. The location 
vector is used as before in conjunction with a dependency vector 
and three tables as shown in Fig. 4. 



Fig. 4 Assembling with Multipoint Constraints 

The dependency vector contains zero entries for all indepen- 
dent displacement components. All dependent displacement compo- 
nents are indicated by a zero in the appropriate location in 

t h 

The index in points to the "d^" location of three 

tables. The first gives the degree of freedom of the first in- 
dependent displacement component in Eq. (39) and the second gives 
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its associated coefficient. With this information, the matrix 
elements in [k] associated with the dependent displacements are 
appropriately multiplied by a^. and then assembled according to 
nj^ in the same manner outlined for single point constraints. 

The last table contains an index that points to the location of 
the a^j and n^ for the next independent pair in Eq. (39). In 
this case, the assembling operation continues as before. A zero 
entry indicates the end of this dependency relation. 

3 . 7 Equation Solver 

Two solution procedures are employed in the PLANS system for 
solution of the matrix equations. They are called PIRATE and 
PODSYM. PIRATE uses a solution algorithm based on the Cholesky 
decomposition scheme and relies on the fact that the matrices are 
positive definite and symmetric, and are explicitly in block tri- 
diagonal form. Because the matrix must be in block tridiagonal 
form, for the most efficient storage allocation, the matrix must 
be partitioned. As a result, all the nodes in each partition 
must be specified and at any one time at least one diagonal and 
one off-diagonal block must be in core. Only the REVBY module 
of PLANS uses PIRATE because of the tight banding possible in 
most axisymmetric shell and body of revolution problems. In ad- 
dition, at most, REVBY considers two dimensional meshes of ele- 
ments which are relatively easy to partition. A detailed de- 
scription of the PIRATE algorithm is presented in this section. 

The second solution procedure, PODSYM, is used in the re- 
maining modules of PLANS. This solution package requires the 
system of matrix equations to be banded, positive definite and 
symmetric. It, too, is based on the Cholesky algorithm but since 
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this algorithm requires only one row of the matrix to be in core 
at any one time, no partitioning is required. Node order must 
be specified, however, and should be ordered so as to minimize 
the bandwidth and thus reduce solution time. A bandwidth opti- 
mizer is available in the SATELLITE program of PLANS for input 
data checking and plotting. 


3.7.1 PODSYM - Solution of Symmetric Positive Definite 
Banded Matrix Equations 

PODSYM solves the matrix equation AX = Y where 

A = 

is a banded positive definite symmetric matrix, X is the de- 
sired solution vector, and Y is the known right side (load 
vector). PODSYM is the user interface and supervisory routine. 

It uses the Cholesky algorithm which factors the total stiffness 
matrix into LL (where L is a lower triangular matrix) and 
then solves a pair of triangular sets of equations. 

The large positive definite matrices that occur in practical 
work very often contain a large number of zero entries and the 
program seeks to benefit from the presence of these elements by 
modifying the standard Cholesky formulas 
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kk 



k-1 



j=l 
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2 


and 
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The factorization is carried out by subroutine QFACT which 
supervises the storage and data set allocation and subroutine 
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QCHOL which generates the lower triangular L matrix. QCHOL im- 
plements the Cholesky algorithm to factor the positive definite 
symmetric A matrix as the product of a lower triangular matrix 
with its transpose 

T 

A = LL 

A straightforward argument establishes the possibility of decom- 
posing any positive definite matrix in this fashion. Once L 
has been obtained, it is not difficult to solve the system of 
linear equations AX = Y by calculating Z as the solution of 

the lower triangular system LZ = Y, and then X as the solu- 

T 

tion of the upper triangular system L X = Z. The forward solu- 
tion (LZ = Y) is accomplished by subroutines QFSOL and QFOR 
and the back solution (L X = Z) by subroutines QBSOL and QBAC . 
Before the call to QBSOL, subroutine REVERS is called, which re- 
verses the rows of L so that the last row becomes the first 

row, etc. This is accomplished in order to sequentially access 
T 

L during the back solution. 

The above algorithm is noteworthy for its assured stability 
and general efficiency. It is possible to carry out an error 
analysis of the procedure as it is represented on a digital com- 
puter. Such analysis shows that the computed L matrix satis- 
fies an equation of the form 

A + E = Ll'^ 

with bounds on the elements of E which show that E is small 

compared to A. The effect of rounding errors made in the sub- 

T 

sequent solution of LZ = Y and L X = Z may then be taken into 
account by (implicitly) introducing an additional perturbation F 
into A, and it is then concluded that the computed solution X^ 
exactly satisfies the equation 
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(A + E + F)X^ = Y 

Since E + F can be shown to be small, one would like to infer 
that is almost an exact solution of the original equations, 

and unless A is too nearly singular, such a conclusion is in- 
deed warranted. But if A is very ill-conditioned, no such re- 
sult can be guaranteed, and X^ may be far from the mathemati- 
cally correct solution; in this event single-precision computa- 
tion will not suffice for the calculation of an accurate solution 
and since the solution will be very sensitive to small errors in 
A, it is unlikely that even a high-precision computation will 
yield satisfactory results unless A and Y are known (and sup- 
plied) to more than single-precision accuracy. The PODSYM sub- 
routines make a fairly realistic attempt to detect and report 
pathological conditions of this sort. 

3,7.2 PIRATE - Solution of Symmetric Positive Definite Matrix 

Equations in Block Tridiagonal Form 

PIRATE solves the matrix equation AX = Y, where 



is the block tridiagonal, positive definite, symmetric stiffness 
matrix and X and Y (representing the generalized nodal dis- 
placements and loads, respectively) are partitioned correspond- 
ingly as 
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This algorithm factors the total stiffness matrix into the pro 
duct of a lower triangular array and its transpose and then 
solves a pair of triangular sets of equations. This factoriza 
tion is possible only for positive definite matrices. 


The method makes use of the Cholesky algorithm to factor the 

matrix A into the product of an upper and lower triangular 

T 

matrix such that A = LL with 
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Ml 4 
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These equations are used in turn to determine L^, L 2 , 

^ 2 , etc., obtaining each diagonal block by the Cholesky algo- 
rithm and each off-diagonal block by solving triangular equa- 
tions . 


The diagonal blocks of the A matrix, A^, are factored 
using the standard Cholesky formulas for each block 


and 


^kk 
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for 


i > k 


A flow chart for PIRATE is shown below. 
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The factorization is carried out by subroutine ATTACK as 
follows 
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^ 2^2 “ ^2 ' 

"" ^3 “ ^ 2^2 


etc . 


MORGAN calls SKULL to form the products and TRIEQ to solve 

the triangular equations for the Zy^. 

The "back" solution is so called because it obtains the 
solution elements in reverse order; in partitioned form the 
process is 

JT Mt 
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2 


X 3 


"3 


^4 


\ ^ J 


1 .e . 


Hh = 2, 


LjXj = Zo - m;x. 




“^3 


etc . 


Once again, SKULL generates the products ^^^+1 TRIEQ pro- 


vides the solutions X, . Since the 

k 


X^ are obtained in reverse 


order, it is necessary for MORGAN to read them backwards in 
order to produce the solution in normal order in core. Note 
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that this is done only after the initial solution since a "back- 
wards" copy of the blocks of L are written on auxiliary stor- 
age for use in subsequent solutions that do not require the 
stiffness matrix to be factored. 


4. ELEMENTS IN THE PLANS LIBRARY 

In this section details about the theoretical aspects of the 
finite elements available in PLANS and their capabilities are 
presented. The order of discussion is based on the analytical 
and/or spatial dimensionality of the elements, i.e., one dimen- 
sional elements are discussed first, then two dimensional ele- 
ments, and, finally, three dimensional elements. The topics dis- 
cussed for each element are presented in the following sequence: 

• Introduction 

• Displacement Assumptions 

• Formation of Stiffness Matrix 

• Geometric Stiffness Matrix 

• Initial Strain Stiffness Matrix (Plastic Load 
Vector) 

• Stress Calculations 

• Thermal Stress Calculations 

• Material Properties 

• Mechanical Loads 

• General Comments 
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4 . 1 Axisypimetric Doubly Curved Shell Element 



Displacement Assumptions: 
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Initial Strain Distribution: 
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Fig. 5 Axisymmetric Thin-Shell Element 
4.1.1 Introduction 

This element is a doubly curved isoparametric thin shell 
element for the analysis of axisymmetrically loaded thin shells 
of revolution. It should be used for shells vjhose thickness to 
mean radius ratio is less than 1/10. The shell theory used to 
derive the element is Sanders' nonlinear theory for small strains 
and moderate rotations (Ref. 4^. The element is based on the 
formulation initially presented by Levine and Armen (Ref. 44), 
and extended in Ref. 26. It has six degrees of freedom at each 
node and an axisymmetric torsion capability has been included. 
Thicknesses are specified at nodes and allowance is made for a 
linear thickness variation from node to node. 
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Tvjo local coordinate systems are employed in the element 
derivation in addition to the global cylindrical r,0,z system. 
One is a local rectangular Cartesian system as shown in 

Fig. 5 with displacements u^, \x^, v. The other is the shell 
curvilinear coordinate system defined by the middle surface dis- 
placement in the meridional direction, u, circumferential di- 
rection, V, and normal direction, w. Here C is the dis- 
tance normal to the middle surface in the direction of w. The 
sign convention for stress and moment resultants and applied 
pressures is illustrated in Fig. 6. 


4.1.2 Displacement Assumptions 

The displacement field may be written in terms of the 
generalized Cartesian nodal displacements as follows 

1 J 1 J 

i J ^ J 


where (O » j (O cubic Hermitian 
interpolation polynomials. This can be written in matrix form 
as 


[u^} = [N][u^J 


(42) 


This results in six degrees of freedom per node. These are 
the Cartesian displacements and their respective first deriva- 
tives. These can be related to the shell displacements, u,v,w, 
the rotation x = dw/ds - u/Rj^, the linear meridional membrane 
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strain = du/ds + and the linear membrane shear strain 

= dv/ds - V cos through the matrix [T] 



( 43 ) 
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or 


{u } = [T]{u 1 
c s 

The final nodal degrees of freedom {u ) are u,w,x>€°,v, 
o s s 

and leading to a 12 degrees of freedom element. These 

quantities are assembled in the tangential-normal coordinate 

system at each node. 


4.1.3 Formation of Stiffness Matrix 

Based on Sanders' theory (Ref, 43) the linear components 
of the strain-displacement relations are 

T ^ + 2Ck „ 

's0 sd s6 

where 


o _ ^ , ]±_ 
^s ~ ds 


e° = — (u cos cp+w sin cp) 



dv V cos 9 
ds r 

o 


K 

S 


dX _ cos 9 

ds ’ ^^0 ~ r 


_ 1 |1 ^ 

s0 ~ 2r^ |2 ds 


/o • Ov , V cos 9 

(3 sin cp - — ) + 


2R, 


3R- sin cp . 

(1 - ) 


X 


dw _ 
ds Rj^ 


The strain displacement relations may now be written in 
terms of the rectilinear displacements u^,U 2 jV and these are 
related to the nodal degrees of freedom to yield 
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( 45 ) 


{€} = (€^3 + aK) = c[w^] + c[w^])(Ug3 = [w]i:ug) 

The elastic stiffness matrix may then be formed as 


[k“] 


[W]^[E][W]dv 


V 


(46) 


= 2w£ 




where [E] is the matrix of elastic constants relating stress to 
strain. The quantity h is the shell thickness which varies 
linearly from node to node. 

This line integral of Eq. (46) is evaluated numerically 
using Gauss-Legendre integration whose order is specified by the 
user. Up to eighth order integration is available but sixth 
order has been found to be sufficient for most purposes. For 
simple geometries and constant thickness, third order is adequate. 

The actual geometry between nodes i and j is approxi- 
mated by a substitute curve, represented in local ^ ,ri coordi- 
nates as Tj = T) (I) . This curve matches the coordinates and 
slopes of the actual curve at the nodes of the element and is a 
cubic Hermit ian polynomial 

q = (2e^- + + (3^^~ 2e^)^2 + (^^- + - 4^)q2 

(47) 

The terms in parenthesis are the first order Hermitian poly- 
nomials Hq^(0» Hq 2^(4)> H^^^(4), H^2^(0* Here and ^2 
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are the ordinates of the substitute curve at the end points and 
r|j^ = T)2-0* The terms and r]^ the slopes tan and 

tan p j . All appropriate geometric quantities such as rj, R, 
and r^ are obtained from the R,z coordinates of points i 
and j and cp^ and cp^ , the angle between the outward normal 
and the positive z-axis. 


4.1.4 Geometric Stiffness Matrix 


The geometric nonlinear stiffness matrix is developed based 
on the small strain moderate rotation theory proposed by Sanders. 
It is intended for use with the "updated" or "convected" coordi- 
nate approach to the solution of geometric nonlinear problems. 

The approximate strain displacement equations specialized to axi- 
S5nnmetric problems are 


% = + Ckg 

•>'se “ •>'s6 + 'Pl”2 + 

where 


u 


dw 

di ^ 


V 


cp = 


1 o ’ 


9 ' 


We can write the geometric stiffness as 


(48) 


where 



[f2]'^[N] [J2]dA 


"A 


(49) 
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( 50 ) 


J cp2>= [nUu^] 

1^3] 


The stress resultants we assumed to vary linearly from node 
to node since stress computations are carried out at nodes 


where 


[N] = [N].(l-0 + [N]^e 


(51) 


This leads to 


[N] = 


N N 


s 6 




0 

0 


0 0 N +N. 

s d 


[k^] = 


[^^]^[N].[n](l-OdA + 


[n] [N]^[n]^dA 


(52) 


where 


dA 


lirr id^ 
o 

cos P 


All area (line) integrations are performed using Gauss- 
Legendre integration of the same order as specified for the 
stiffness matrix. 


4.1.5 Initial Strain Stiffness Matrix (Plastic Load Vector) 
The effective plastic load vector can be written as 
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The initial (plastic) strains are assumed to vary linearly 
from node to node while at the nodes the variation of the plastic 
strains through the thickness is arbitrary. We have then 

[AeP(e,C)]= fA€P(0,O)(l-O + {AeP(l,C)]l (54) 

Using this expression in Eq. (53) and remembering that [W] = 

[Wm] + final expression for {AQ^} is 





ph/2 

(AQ.} = 

(l-4)[W^]'^[E]dA 


[A€P(0,C)}dC 


A 

t 

-h/2 


,h/2 


+ I i[Vl] [E]dA 
m 


(A€P(l,0}d? 


-h/2 



•> 


ph/2 

+ 

(l-0[W^]'^[E]dA 


aAeP(0,O)dC 

f. 

A 

t 

-h/2 





ph/2 

+ 

^[W^]'^[E]dA 


aAeP(l,C)}dC 

t 

A 

t 

-h/2 


(55) 


The integration of the effective plastic forces and moments 
through the thickness is carried out using Simpson's rule with an 


option of choosing up to 21 integration points (20 layers) . 

The number of layers used must be an even number and is input by 
the user. Eleven points (10 layers) are usually sufficient. 

The area (line) integration is performed using a Gauss -Legendre 
integration scheme (input by the user) of the order used to form 
the stiffness matrix. 

4.1.6 Stress Calculations 

Stresses are calculated at the nodes of the element. All 

strain components are continuous from element to element except 

K . This quantity is averaged at nodes. All stresses output 
s 

are in the shell curvilinear coordinate system. Stresses are 
also calculated and output at each integration point through the 
thickness by node. 

4.1.7 Thermal Stress Calculations 

Orthotropic thermal stress calculations are allowed. Dif- 
ferent thermal coefficients of expansion in the meridional and 
circumferential directions can be input. A parabolic temperature 
distribution through the thickness is assumed at nodes. Between 
nodes the meridional temperature variation is assumed to be 
linear. Temperatures are input at the top, bottom, and middle 
surface at each node. The thermal load vector is obtained from 
Eq. (55) where is replaced by a^i^AT, i.e., 

OgATl 

{As^} = (56) 

The integration of the thermal forces and moments is done 
exactly, assuming the parabolic temperature distribution. 
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Temperatures at layers are obtained by interpolating their values 
based on the parabolic distribution. 


4.1.8 Material Properties 


All material properties are constant within the element and 
temperature independent. Elastic properties are orthotropic. 

The meridional and circumferential directions are the principal 
directions of orthotropy. We have 


> 


a 

s 


“e 

II 




1 - V 

s6 6s 


V E 

6s e 
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s6 s 


1 - V V 

s6 6s 


1 - V V 

se es 
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'se 
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S 
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/ 


0 




[C][e] 


where ” ^12 ^ ^ positive definite stiffness matrices 


Plastic material properties are assumed to be orthotropic 
with Hill's yield criterion for plane stress used for initial 
yielding 

f = (G+H)ag + (F+H)0g - 2H0g0g + 2NTg^ = 1 (58) 

with G+H=l/X^, F+H=l/Y^, 2H = 1/X^ + 1/Y^ - 1/Z^ , 2N=1/T^, and 
X,Y,Z are the yield stresses in tension in the s,6,^ direc- 
tions, respectively. T is the yield stress in shear in the s-0 
plane. For two dimensions the stability criterion for use of 
this yield condition reduces to 


i i . 1 (A+ j. 

2 2 4 \ 2 2 

X Y X Y 



(59) 
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Perfectly plastic j linear strain hardening or nonlinear 
strain hardening laws may be chosen. For orthotropic plasticity 
the individual components of the stress-strain law are input and 
each component may be represented as linear or nonlinear strain 
hardening. If one component is perfectly plastic all must be. 

At the end of each half-cycle of loading new plastic material 
properties may be input but not new elastic material properties. 
This allows subsequent half-cycles of the stress-strain curve 
to be accurately represented. 

Thermal coefficients of expansion may be input orthotropi- 
cally, i.e., different in the meridional and circumferential di- 
rections. They may not be changed at the end of each half-cycle. 

4.1.9 Mechanical Loads 

Two types of mechanical loads may be applied to a shell ele- 
ment . They are surface tractions and concentrated (or line) 
loads . 

Concentrated or Line Loads (Edge Loads) — These are applied 

as forces per unit length of circumference except at r = 0 

where N =M =N„=0 and F = Q represents an actual 
s s s6 z s 

force. At any other radial location the line (or edge) loads are 

applied in the natural (shell) curvilinear coordinates as N , 

s 

Q , M , and N „ in the u,w,x, and v directions, respec- 

o o o O 

tively. These forces per unit length are applied and specified 
at nodes of the shell. 

Surface Loads — The surface tractions are specified per unit 
surface area in the meridional p , normal, p„ , and circumfer- 
ential directions, p^. They are assumed to vary linearly from 
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node to node. All three components at each node must be speci- 
fied. The loads are applied through consistent load vectors 
which are calculated in the program 




P 


1 

fPs. - Ps.ll 




X 




[P^] = 2ir£ 


[Nj'^rdl^ 


>+ 

[Ni'^erde^ 

P^. ■ Pc. ) 







J 1 



0 


t 

0 







J 1^1 


All integrations are carried out numerically using a Gauss - 
Legendre scheme of the same order as specified for the stiffness 
matrix. 

4.1.10 Equilibrium Correction 

The equilibrium correction term is written as 


= [P^} - 


[W] [ojdv = {P. } - 


[W]^(l-0(N^}dA 


(61) 


[W]^^[Nj}dA 


Here the stresses are assumed to vary linearly from node to node 
for the equilibrium correction. The stress and moment resultant 
vectors at each node {N^} are evaluated by numerically inte- 
grating the stresses through the thickness using Simpson's rule 
in the same manner used to obtain the "inelastic forces and 
moments ." 
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4.1.11 General Comments 

• Pole Conditions — At a pole (r = 0) the following condi- 
tions must be satisfied for the shell element: 1) = 0, if 

we have a shell where dr/dz = 0 at r = 0 (i.e,, cp = 0 or 
ir) then this implies u = 0; 2)x = 0j 3)v=0; and 

4) = 0. Correspondingly, no forces can be applied at 

r . 0, I.e., 0. 

• Cap Element — No special cap element was required for the 
axisymmetric shell element described here. The use of the Hermi- 
tian polynomial representation of the displacement components 
and U 2 eliminates the necessity of obtaining a relationship be- 
tween the generalized coordinates and the degrees of freedom any- 
where in the element, and specifically at a pole. Furthermore, 
since the integration was performed numerically using the Gauss- 
Legendre technique, no actual evaluation of the function matrix 
[W] at the point of singularity or indeterminacy was required. 

The singularity was isolated and all improper integrals required 
for the various matrices converge. The terms of the matrix [W] 
required for the u^ and degrees of freedom, which might 

result in numerical problems, were removed since the radial dis- 
placement and rotation are identically zero at the pole. To ob- 
tain strains at the pole, appropriate limits were taken in the 

[W] matrices, employing L'Hopital's rule. It is, however, neces- 
sary to use a finer grid in the region of the pole than elsewhere. 
Treating the pole as if it were an actual physical boundary leads 
to accurate results for displacements and stresses. 

• Geometric nonlinear analysis using this shell element is 
only available in a special purpose module called AXSHEL. It is 
not available in REVBY at the release date of this report. 
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This element is available in the REVBY module of PLANS. 




4.2 Axisymmetric Ring Element 

4.2.1 Introduction 

Classical beam theory forms the basis for the development 
of the stiffness matrix for a thin circular ring of arbitrary 
cross section, as derived in this section. This element is used 
in the REVBY program and is intended for use with the revolved 
axis 5 nnmetric shell element. Pertinent geometry defining the 
ring-shell intersection is shown in Fig. 7. In the following 
derivations it is assumed that the shear center and centroid of 
the ring coincide. 

4.2.2 Displacement Assumptions and Formation of Stiffness 
Matrix 

Based on thin ring theory for an axisymmetric circular ring, 
the total strain and stress for any point in the cross section in 
the presence of an initial (plastic) strain is 

e = " ~ (u^ - yP^) ( 62 ) 

c 

a = E(e - e) (63) 

where r^ is the radius of the ring, e is an initial (plastic) 

strain, u^ is a component of displacement of the ring axis, and 

P is a rotation normal to the ring cross section, 
z 

Substituting Eqs. (62) and (63) into the expression for 
strain energy leads to 
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( 66 ) 


EA i 0 
[k] = ~ _ - ' _ - 

^0 .El 

* X 

and [k ] is the initial strain stiffness matrix. 

4.2.3 Initial Strain Stiffness Matrix 

The initial strain stiffness matrix obtained in Eq. (65) is 
given by 

[-1 I 0 

[k*] * 2irE — (67) 
Oil. 

The plastic load vector {Q} is the product of [k ] and the 
vector of integrated plastic strains. 

Since the plastic strain irfay vary through the ring cross 
section the integrals in Eq. (64) cannot be determined a priori. 
These integrals are performed numerically for plastic rings using 
a Gaussian quadrature. 

4.2.4 Transformation to Shell Coordinate System 

The stiffness matrices of Eqs. (66) and (67) are written 
with respect to the ring coordinate system. In order to use the 
ring element with the axisymmetric shell element the ring must 
be related to the shell coordinate system. 

The pertinent geometry defining the ring-shell intersection 
is shown in Fig. 8 for an attachment at point A. Continuity of 
the ring and shell displacements at point A leads to 
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where 


( 68 ) 


r 

-r 


[T] 





2 r' 
o 


I 

0 I -1 


^ *\f <\j 

u, w, <t>^ are the shell displacements and rotation, r and R 2 
are the shell radii of curvature shown in Fig. 8, and "prime" 
denotes differentiation with respect to the arc length, s. 



Fig. 8 Geometry Defining the Ring-Shell Intersection 
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The final stiffness and initial strain matrix with respect 
to the shell displacements then is 


[kg] ' [Il’^tkHT] 

(69) 

[k*] = [T]’^[k*l 

(70) 


4.2.5 Stress Calculations 

Stresses are monitored at a specified number of Gauss points 
for each cross section via Eq. (63). From these stress values, 
each point is checked for yielding and if necessary plastic 
strains are calculated. Once these values are determined the 
numerical integration can be performed. 

Five distinct ring cross sections are currently in the REVBY 
module, a solid rectangle, solid circular, Z-section, I-section, 
and hollow circular section (see general comments for section de- 
tails) . 


4.2.6 Thermal Stress Analysis 

A parabolic temperature distribution can be specified through 
the depth of the ring (i.e., along x in Fig. 8). Three tempera- 
tures are specified. They are at the top, centroid, and bottom 
of the ring. The thermal load vector is formed in the same manner 
as the plastic load vector with e replaced by aAT. Thermal 
strains at integration points are determined by interpolation 
based on the assumed parabolic distribution. 

4.2.7 Material Properties 

Elastic, ideally plastic, linear strain hardening or non- 
linear strain hardening stress-strain laws may be chosen. The 
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nonlinear hardening is based on a Ramberg-Osgood representation 
of the actual stress-strain data. The plastic material proper- 
ties may be changed at the end of each half-cycle of loading. 

4.2.8 Loads 

Only axisymmetric line loads may be applied at the node de- 
fining the ring location. 


4.2,9 General Comments 

• Five different cross sectional shapes for the ring are 
available in the REVBY module. These are shown in Fig. 9. They 
are: 1) a solid rectangular cross section, SREC; 2) a solid 

circular cross section, SCIR; 3) a Z-section with equal area 
flanges, i.e., a = b, t^ = t 2 > ZSEC; 4) an I-section with 
equal area flanges, i.e., a = b, ~ ^ 2 ^ ISEC; and 5) a hol- 
low circular cross section, HCIR. The order of the Gauss- 
Legendre integration schemes used to obtain the plastic load 
vector follows 


Section 

SREC 

SCIR 

ZSEC 

ISEC 

HCIR 


Gauss -Legendre Integration Order 

Three by three array (9 points) 

Two radial points times four circumferential 
points (8 points) 

Third order in each flange and web (9 points) 

Third order in each flange and web (9 points) 

Eighth order around the circumference 
(8 points) 
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• For the Z- and I-sections the thicknesses are assumed to 
be small compared to the flange and web lengths. For the hollow 
circular section it is assumed t/a « 1, 

• The applicable module for this element is REVBY. 


4.3 Stringer Element 
4.3.1 Introduction 

This element is used to represent a one dimensional axial 
force structural member. Two stringer elements are included in 
the PLANS system: a two-node element developed from a constant 
axial strain assumption and a three-node element developed from 
a linearly varying strain assumption (Fig. 10). 


4.3.2 Displacement Assumption 


The distribution of the three local displacements 
shown in Fig. 10 is assumed to be 


'’u(x)^ 


u. 

1 


r ^ 

u. 

J 

( V(x) 

>= (1-0^ 

V. 

1 

>+ 


w(x)^ 


L 1 J 


w . 

J -7 


u , v,w 


(71) 


where | = x/i for the two-node stringer, and 


u(x)" 


fuO 


J 



<v(x) 

>= (24^-3^+!)^ 

V. 

X 

y+ (21^ -?) / 

V. 
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>+ 4(?- 
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V7. 
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w. 
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for the three -node stringer. 
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(a) Constant strain element (b) Linear strain element 

Displacement Assumption: 

j 2 

u = a^ + a^x u - a^ + 0 . 2 ^ i a^x 

Initial Strain Distribution 
e = constant 

Fig. 10 Stringer Element 

As shown in Fig, 10, the indices i,j denote the end point 
nodes of the three-node stringer. The specification of the axial 
displacement component is sufficient to fully describe the elastic 
linear response of a stringer element since the element has no 
out-of-plane (normal to the axial direction) stiffness. The nor- 
mal components, v,w, contribute to the elastic, geometric non- 
linear response by coupling the axial force to the out-of-plane 
response. This coupling leads to the initial stress stiffness 
matrix [k^] . It should be noted that displacement components 
v,w, as shown in Fig, 10, refer to any two mutually perpendicu- 
lar directions normal to the stringer axis since the stringer 
does not have a preferential cross sectional reference system. 
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4.3.3 


Element Stiffness Matrix 


The element stiffness matrices are obtained from the princi- 
ple of virtual work (Section 3) . The appropriate strain displace- 
ment relations are 


_ ^ 





J. 


(73) 


Using this approach, the linear and nonlinear components of total 
strain are obtained from Eqs. (71) and (73) as 


Ae = i [-1 I 1] 
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(74) 


(75) 


for the two-node stringer and 


Ae = i [4^-3 j 4?-l I se-4] 


u. 
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“J 


u. 


(76) 
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for the three-node stringer. 

Making use of Eqs. (74) and (76) in the definitions for the 
elastic stiffness matrix [k ] and assuming a constant cross 
section stringer then yields 

1 -1 


tk°i = f 


-1 


for the two-node 
stringer 


(78) 


and 


[k°) = 


7 

-8 

1 


-8 

16 

-8 


1 

-8 

7 


for the three -node 
stringer 


(79) 


where A is the section area and E is Young's modulus. 


4.3.4 


Geometric Stiffness Matrix 


Using Eqs. (75) and (77) in the definition of the geometric 
stiffness [k^] gives 


[kl] = E 

t 
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-1 
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-1 
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0 
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-1 


0 
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-1 

1 


for the two-node 
stringer 


(80) 
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and 


■ 7 1 -8 0 0 0 

1 7 -8 0 0 0 

~ -8 -8 16 0 0 0 

[k^] = J 7 Q Q Q 1 I .Q for the three-node (81) 

stringer 

0 0 0 1 7 -8 

.0 0 0 -8 -8 16 

Here F is the average axial force. 

4.3.5 Initial Strain Stiffness Matrix 

The initial strain matrix is based on the assumption that 
the thermal (plastic) strain is constant in the element. Thus, 
the initial strain matrix can be written 

- 1 ' 

[k*] = EA +1 (82) 

. 0 , 

for both the two- and three-node stringer elements. 

4.3.6 Stress and Force Calculations 

Element axial forces in the presence of thermal (plastic) 
strains are calculated as 

ff) = [k°] f6) - [k*lf£p) (83) 

where {6} are the local axial displacements and {Sp} are the 
average thermal and/or plastic strains. Average element stresses 
are calculated at the center of the element and are used in conjunc- 
tion with plasticity calculations consistent with the assumption of 
constant plastic strain within each element. 
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4.3.7 


Thermal Stress Calculations 


Temperatures are specified at the two end point nodes. 

These nodal values are used to determine an average element 
temperature . 

4.3.8 Material Properties 

Plastic strains are calculated using a uniaxial bilinear 
stress-strain relation or a nonlinear Ramberg -Osgood relation. 
Perfect plasticity is also accommodated. 

4.3.9 Loads 

The stringer can be loaded by concentrated nodal loads or 
consistent distributed line loads (Section 4.7.9). 

4.3.10 Comments 

• The three-node stringer has an assumed displacement function 
consistent with that used for the LST of Section 4.7 and is, 
therefore, used when linking an LST with a stringer. 

• The applicable modules for this element are OUT-OF-PLANE, 
BEND, and OUT -OF -PLANE -MG. 

4.4 Beam Element 
4.4.1 Introduction 

The stiffness properties for an initially straight beam 
finite element of arbitrary cross section are derived in this 
section based on the assumptions of classical beam theory. These 
assumptions include 
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1 ) 


Normals to the centroidal axis remain straight and 
normal after deformation and their length remains un- 
changed, i.e,, the effect of transverse shear deforma- 
tion and transverse normal strains are neglected. 

2) Warping of the cross section is nelgected. 

Although assumption 2) is restrictive for all but circular cross 
sections, the ability to specify the actual torsional rigidity 
has been maintained. Throughout the development, it is assumed 
that the shear center and centroid of the cross section do not 
necessarily coincide. The element coordinate system (located at 
the shear center of the cross section) and convention for mo- 
ments, forces, and displacements are shown in Fig. 11. Based on 




Fig. 11 Coordinate System and Convention for Forces 
and Displacements for Three Dimensional 
Beam Element 
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the convention of Fig. 11 and assumptions 1) and 2), as noted 
previously, the displacements at any point within the beam cross 
section can be written as 


= u + (z - 25 )P ” (y " y )P 
= Vs - (84) 


w = w + yP 
z s X 

where u^ is the axial displacement of the centroidal axis; v^, 

w^^ are the lateral displacements of the shear center in the 

cross sectional y and z directions, p , p , p are cross 

y Z 

sectional rotations defined in Fig. 12, and y ,z are the 

® -^cg eg 

distances between the centroid and shear center. 




z 



Fig. 12 Convention for Rotations 
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Based on Eqs . (84) and excluding terms associated with non- 
linear terms in the curvatures, the nonlinear strain displace- 
ment relations are 


e = e + (z - z )k " (y “ y )k. - 

X o'* cg‘^ y -^cg' z X 


7 =-66 - ZK 


xy 


X' y 


-Tf 


yz ’ xy 


(85) 


y = yK - p p - 7^ 


yz 


X z xz 


where 


■ SX + 2^2 


K 


Sp 


X 


yz ^x 


K 


" ■ Sx2 


and e^, 7 ^ , 7 ^ are initial or plastic strains. 

X xz ’ ' xy ^ 

Equations (85) along with the definitions for moment and 

product of inertia about the section centroid 


I = 

y 


(z - z ) dA , I = 
^ cg^ ’ z 


(y - y<;g) dA 


A 


yz 


(y - yeg)(2 - ^cg^*** 


are used to derive the element stiffness matrices 
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4 . 4.2 


Displacement Assumptions 


There are four independent displacement components that are 

prescribed to formulate the beam element stiffness properties. 

These are the two transverse displacement components, v , w 

s s 

and the axial displacement and twist, u and p , The out-of- 

O 

plane rotations are related to the transverse displacements in 
the usual manner by imposing Kirchoffs hypothesis, i.e.. 


P„ = -■« 


X 


= V 


X 


the two components of transverse displacement are assumed to be 
cubic functions that are solutions to the linear homogeneous beam 
equations 






V 


r I ' + 

T = 1 I V7 . * 

^ 1 




w i=l 

s ' 


I w 

I .’=1 

where ^ = x/i, i is the element length, and 

- 31^+1 

= -2e^ + 3^^ 


( 86 ) 


It is assumed that the in-plane displacement and twist u , P 

C X 

are linear functions of the axial coordinate 
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(87) 


where 


(“cl 


2 


- I 



i=l 


i 




H 


( 0 ) _ . 

02 ^ 


4.4.3 Stiffness Matrix 


The elastic component of the stiffness matrix can be de- 
veloped using the principle of virtual work or from energy con- 
siderations. In either case, use is made of the kinematic dis- 
placement assumptions, Eqs . (84), to reduce the general relation- 
ship for the stiffness matrix, which is a function of the three 
spatial coordinates and involves a volume integral, to a function 
of only the axial coordinate and a simple one dimensional inte- 
gral of the form 



J 

[W]'^[D] [W]dx 


^ 0 


where 



EA 

0 

0 

0 ■ 


0 

GJ 

0 

0 

[D] = 

0 

0 

El 

yy 

-El 

yz 


0 

0 

-El ^ 
yz 

El 

zz 1 
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and A is the cross sectional area, I , I , I the cen- 

yy’ zz’ yz 

troidal moments and product of inertia, respectively, J is 
the section torsional rigidity, and E and G are Young’s 
modulus and the shear modulus, respectively. The [W] matrix 
is obtained by substituting the displacement assumptions, 

Eqs. ( 86 ) and (87) into the linear component of the strain dis- 
placement relations, Eq. (85) and is written as 

f--jO 0 0 00^ oooool 


0 0 0 


1 


000 OO 7 OO 

£ 


[W] = 


00 --“ 0-—0 0 0 ^0 


0 


( 88 ) 


<}> ^2 '“’l ''’3 

0^0 0 0^0-^000-^ 

where 

= 6(24 - 1) , <t>2 = 

<1>3 = (64 - 2) 

The resulting stiffness matrix is a 12 x 12 matrix with 

six degrees of freedom per node u.,v.,w.,p ,3 , 

111 x^ y^ z^ 

where the axial displacement refers to the centroid and the 
lateral displacements refer to the shear center. 


4.4.4 Geometric Stiffness Matrix 

The development of the "geometric" stiffness matrix or ini- 
tial stress stiffness matrix, based on small strain, moderate ro- 
tation assumptions for a planar beam, is straightforward and has 
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r 


been previously presented in Ref. 17. However, for the case of 
a general space frame subjected to torsion as well as bending, 
there is nonlinear coupling between bending and torsion which can 
lead to various expressions for the initial stress stiffness 
matrix depending on the approximations made in arriving at the 
strain displacement relations, Eq. (85). In this development 
the simplest approach was taken in that nonlinear contributions 
to the curvature and twist are neglected and Kirchofffe hypothesis 
is imposed on only the linear components of strain. In so doing 
the form of [k^] is similar to that for a planar beam with 
the addition of only the nonlinear coupling between torsion and 
bending. 

With these considerations, [k ] can be written as 


where 


and 




[k^] = 


[af [2][n]dx 


(89) 


[ 2 ] = 


0 -V -V 


-V 


-V 


T 

0 


0 

T 


0 0 0 < t >, 000 00 ( 1)^00 


^1 

[fi] = JO 0 ’ I ^ < 1>2 0 0 0 <t>j^ 0 <*>30 


Ot^ 0 0 0<t>^0-~0 0 0<|)v 
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4.4.5 


4.^ = - O , 4>2 = 34^ - 4^ + 1 , 4)3 = 34^ - 24 

0>4 = 1 - I , ^5 = ^ 


Initial Strain Stiffness Matrix 


The initial strain stiffness matrix is used to compute the 
nodal forces due to thermal or initial strains. In a plastic 
analysis using the "initial strain" method it is used to deter- 
mine the effective plastic load vector (Ref. 8) as 

[AQ] = [k*HAe} 


where 

{AQ} is the effective plastic load vector 
[k ] is the initial strain stiffness matrix 

and 

{Ae} is the vector of discrete quantities character- 
izing the state of plastic strain. 


The initial strain stiffness matrix, based on the displacement 
assumptions of Eq. (84) can be written as 


[k*] = 


Ae 


X 




(z - z )Ae 

^ cg^ X 


(y - ycg)'^"x 


) dV 


P P 

vAy^ - zAy 
^ ^xz ^ yz 


(90) 


where [W] is the matrix defined in Eq. (88) and 
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0 


0 


0 


FE 


[E] 


0 E 0 0 

0 0 -E 0 


0 0 0 G 


Because for a three dimensional space frame, the plastic strain 
varies in each cross section as well as in the axial direction, 
Eq. (90) cannot be immediately reduced to a one dimensional in- 
tegral as was the case for the stiffness matrix. However, con- 
sistent with finite element analysis, the functional form for 
the distribution of plastic strain increment can be assumed for 
each beam element. To this end, it is assumed that each compo- 
nent of increment of plastic strain varies linearly in the beam 
axial direction. Thus 


"^xy = 

1 J 

where the i,j denote the two beam nodes, respectively (Fig. 11). 
At this point no assumption is made for the distribution in any 
cross section. With this assumption, [k J can be written as 
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[k*] - - 


[W]'^[E][C, 


where 


and 


0 


C2ldx 


[C^] = 


i = 1, 2 


X 


0 


0 


0 


Ae 


X, 


Ae 


X. 


)dA 


P P 

yAy^ - zAy 

xz. xy.^ 


0 


X 


0 


0 


0 


0 


X 


0 


0 


0 


0 


1 


(92) 


4,(0) = (1 . ^) . <t,(<^) = ^ 


The integrals in Eq. (92) are evaluated numerically using a 
Gauss -Legendre integration scheme. To accomplish this, the 
shape of the cross section must be known a priori and increments 
of plastic strain must be evaluated at the Gauss points in the 
cross section. To this end nine distinct cross sections have 
been provided for in PLANS. These are shown in Fig. 13 along 
with the order of Gauss -Legendre integration. 
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Solid Rectangular Section 
l6 Integration Points 


Solid Circular Section 
8 Integration Points 



Z ~ Section 
9 Integration Points 



Hollow Circular Section 
8 Integration Points 



L - Section 
6 Integration Points 



I - Section 
9 Integration Points 



Hollow Rectangular Section 
12 Integration Points 


^2 



T - Section 
6 Integration Points 



y 


Channel Section 
9 Integration Points 


Fig. 13 Beam Sections Available in PLANS 
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4.4.6 Transformation to Global Coordinates and Offset 
Specification 

Once the element stiffness matrices have been evaluated, it 
is necessary to assemble each matrix to some common global sys- 
tem. However, since the beam element is geometrically one dimen- 
sional, only the orientation of the beam axis is known by speci- 
fying the coordinates of the end point nodes. Consequently, the 
orientation of the cross sectional axes remains undefined. In 
order to define this orientation an extra node point must be 
specified as shown in Fig. 14. 

This additional node defines a plane that fixes the direc- 
tion and orientation of the y-axis. The z-axis is perpendicu- 
lar to this plane. In practice, this node can be part of the 
structure or it can be specified only for the purpose of defining 



Fig. 14 Beam Node Specification 
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the cross sectional orientation. If the latter option is chosen 
all degrees of freedom associated with the node must be fixed. 
Once this additional point is defined, the direction cosine 
transformations between the local and global directions are com- 
pletely determined. 

Also accounted for in the transformation between the local 
and global systems is the offset of the attachment point (point 
A in Fig. 15) from the beam axis. The appropriate transforma- 
tion is obtained from the beam kinematic assumptions, Eq. (84), 

t tl 

and can be written for the i node as 



Fig. 15 Definition of Offset Point 
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Geometric nonlinearity is considered in PLANS using an up- 
dated coordinate system. In using this method the change in 
orientation of the beam element must be accounted for during 
each load step. To this end the change in orientation of each 
element is calculated as 









m - 

a 6 










-<^e”a + 
■*a”e“p) 



where the coordinates shown are with respect to the beam local 
axis system and 
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ig = cos b 

= cos ^ 

P 

i = cos a 

a 


nig = sin 6 

= sin p 

P 

m = sin a 

a 


Superscripts i+1, i represent the current and previous con- 
figuration, respectively, and £ is the element length in the 
i^ configuration. 


4.4.7 Force, Moment, and Stress Calculations 

The incremental element nodal forces and moments are ob- 
tained in each load increment by making use of the element stiff- 
ness matrices related to the local element coordinate system. 

Thus 

{Af3 = j[k^] + [k^]|{A6^3 - [k*UAe3 (93) 

Total forces are obtained by summing incremental quantities. 

The stresses are evaluated within each increment by substi- 
tuting the displacement assumptions, Eqs . (84), into the linear 
portion of the strain displacement relations and then making use 
of the appropriate stress-strain relations. These equations are 
evaluated at the Gauss points of each section. 
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4.4.8 Material Properties 

Elastic, ideally plastic, linear strain hardening, or non- 
linear strain hardening stress-strain laws may be used with the 
beam element. The nonlinear hardening is based on a Ramberg- 
Osgood representation of the actual stress-strain data. The 
plastic material properties may be changed at the end of each 
half-cycle of loading. 

4.4.9 Loads 

The beam element can be loaded by 

a) Concentrated forces and moments at each node 

b) Linearly varying distributed load in the local 
y and z directions. 

The consistent load vector for the distributed load is obtained 
from 


where [N] is the shape function for the lateral displacement 

- 3 ^^ + 1 
- 2 ^^ + 3 ^^ 

[N] = (94) 

3 2 

i(r - r) ) 

[Q] = [1 - e ! e] 
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are the magnitude of the distributed load at each node 
in the y or z direction and Cf} is the vector of local forces 
and moments in the x-y or x-z plane. 

4.4.10 Comments 

• The beam element is available in the OUT-OF-PLANE, OUT-OF- 
PLANE -MG, and BEND modules of PLANS 

• Currently, offsets at both ends of the beam must be the 
same 

• No thermal capability is currently available for the beam 
element 

4.5 Shear Panel Element 

4.5.1 Introduction 

A warped shear panel has been included in the element 
library of PLANS. This element is used in practice to model 
shear webs of wing spars and ribs as well as the sheet area of 
some fuselage sections. The development of the stiffness proper- 
ties for this element does not follow the usual consistent de- 
velopment of element properties. The procedure shown herein fol- 
lows the development as described in Ref. 46. 

4.5.2 Stiffness Matrix 

The element characteristics for the warped shear panel are 
developed from the assumption of a distribution of stresses 
within the element that satisfy equilibrium but do not satisfy 
strain compatibility, except in the case of a parallelogram. A 
flexibility coefficient is then computed, using an energy formu- 
lation, and a set of equilibrium equations is used to obtain the 
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stiffness matrix. Details involved in obtaining the energy ex- 
pression can be found in Ref. 47. In order to obtain the stiff- 
ness of the element, six equilibrium equations are written for 
the element in the reference coordinate system shown in Fig. 16. 

The forces f^, f', f^, and f^ (Fig. 17) are necessary to 

ensure equilibrium in the z-direction. The next step is to 
^ ^ ^ ^ 

solve for f^^, f^, f^, q 2 > <^3 > q^, in terms of q^. 

Since there are six equations and seven unknowns, an additional 
equation is needed, which is obtained by assuming that the re- 
sultant force in the z-direction passes through point v 
(Fig. 18), and that one-half of this resultant is acting at 
nodes j and k. 



Fig. 16 Projected Quadrilateral 
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Taking moments about a line passing through node j and 
parallel to i-^ 


g/ 


^3 * i + ^2 + ^3 + 


or 


where 


f' + f' + (1 - 7)f' + f' = 0 


(95) 


(96) 


7 = 


2(83 + g^) 


(97) 


The equilibrium equations and Eq. (96) can be written in the 
matrix form 


[EQ][£”} = [R}q^ (98) 

where [f"} = £3 <?2 ^3 94 ^* Solving Eq. (98) we have 

(f") - [EQ]"’-{R)q^ (99) 


or 


{f} = {E]q^ (100) 

where {f} is obtained by enlarging [f’} to include q, and 
[E] is obtained by correspondingly enlarging [EQ] [R] to in- 
clude the identity q^ »= q^^. 

The eight applied forces can be considered to correspond to 
eight degrees of freedom. They are related, as we have seen, by 
seven equations, six of which are based on static equilibrium 
and the seventh is an assumed relationship. An additional rela- 
tionship, based on elastic deformation and corresponding to a 
single elastic degree of freedom can be introduced. If we regard 
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as representing the generalized force corresponding to the 
elastic degree of freedom and 6 is defined as the correspond- 
ing generalized displacement, the elastic strain energy may be 
written in the form 


U = 


1 . 

2 


Introducing the relationship 


6 q, 
qi 1 


( 101 ) 


6 


q 


1 


aq^ 


( 102 ) 


•where a is a flexibility coefficient, we can write Eq. (101) in 
the form 


U = 


2a 



(103) 


An alternative form for the strain energy is 

U = i {f}'^{6} (104) 

where {5} is a displacement matrix corresponding to {f}. Sub- 
stitution of Eq. (100) into Eq. (104) yields 


U = i [E}'^[6}q3_ 

Comparing Eqs. (101) and (105), we see that 

'll 

Substitution of Eq. (106) into Eq. (103) yields the form 

U = ^ {6}^{E}{E}'^C6) 
or 

U = 4 (6)’’[k]t8) 

where 


(105) 


(106) 


(107) 


(108) 
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I 



(109) 


[k] - i [E)(E)’' 

[k] is seen to be the stiffness matrix in the relationship 

{f} = [k]{6] (110) 

The flexibility coefficient a is computed for the pro- 
jected geometry of the warped quadrilateral on the reference 
plane defined by the lines 1-4 and 2-3 in Fig. 16. This coeffi- 
cient gives the relationship between the generalized shear de- 
formation, 6 , and the shear flow, q, , acting along side 

4l 

1-2, as shown in Fig. 16. 

The computation of a for a particular quadrilateral de- 
pends on its shape, which can be any c£ the following four types 

• Parallelogram or rectangle 

• Trapezoid with side a parallel to side c 

• Trapezoid with side b parallel to side d 

• General quadrilateral 

In general, the expression for a is not easily obtained 
for trapezoidal and quadrilateral shapes. Approximate expres- 
sions have been derived, however, for these shapes by Garvey 
(Ref. 47) based on the assumption that there are some directions 
in the panel that are in states of pure shear. This is true for 
both rectangles and parallelograms, but is only an approximation, 
for trapezoids and quadrilaterals. Therefore, shapes which vary 
substantially from a parallelogram introduce errors in the analy- 
sis . 
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4.5.3 


Initial Strain Stiffness Matrix 


In the presence of initial strains an initial generalized 
displacement component, 6° , must be added to Eqs . (102) through 
(103) so that 


= aq. + 

qj. 1 q, 


and 


6^ = {E}'^([6} - [6)°) 


where (5}° is an initial displacement associated with the 
forces {f} of Eq. (104). This leads to the expression for 
strain energy 

U = i[6}^[k][6] -i[5l'^[E][E]T{6]° + i(6°}T[E][E]^[6°] 


Since 


and 


(6° ) = CE)’^t6°) 

(6°^) = (Ef a 


where is an initial strain corresponding to q^; 

plete integral can be written as 


the com- 


U= i(6}^[k](6}^ - {6}^[k]\ - 4[5°}^[E}(E}^[6°} 
where the initial strain matrix is 

[k ] = [E]^ G 
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4.5.4 


Comments 


• The shear panel is available in the OUT-OF-PLANE module 

• The element can currently be elastic only 

• The element has no thermal capability 

• Force output is in terms of the shear flows and normal 

forces 


4 . 6 Revolved Triangular Ring Element 
4.6.1 Introduction 

The revolved ring element is used to describe three dimen- 
sional axisymmetric bodies subjected to axisymmetric loads. The 
element is two dimensional in that there is no variation of 
stress or strain in the circumferential direction. However, the 
strains and stresses in the circumferential direction must be 
considered as they are induced by radial displacements and 
anisotropy. 

The element is an orthotropic triangular ring element simi- 
lar to that described in the literature (see, for example. 

Ref. 45). Figure 19 depicts the element along with the principal 
material directions. 


4.6.2 


Displacement Assumptions 


The element is characterized by a linear displacement as- 
sumption, independent of the circumferential direction, i.e.. 


u = a, + a„r + a^z 
r 1 2 3 


u = a, + a^r + a,z 
z 4 5 6 




( 111 ) 
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z 



r 


Fig. 19 Revolved Triangular Ring Element and 
Principal Material Directions 


The appropriate strain displacement equations are then 


^rr 

= u 

r,r 

^6Z 

= 

0,z 



= u /r 

G 

= u 

+ u 

^ee 

x' 

rz 

r ,z 


"zz 

= u 

z , z 

CD 

= 

e,r 

- u 


Substituting Eq. (Ill) into Eq. (112) we arrive at 

[e] = [B](a] 


where 

T 

“ [^60 ^zz ^rr ^rz ^r0 ^z0j 

. T 

[a] - ^2 • • * ^9 J 


( 112 ) 


(113) 
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and 


1/r 1 z/r 0 0 0 0 0 0 

0 000010 00 

0 100000 00 

0 010100 00 

0 0 0 0 0 0 -1/r 0 -z/r 

0 000000 01 

From Eq. (Ill) we can determine {a} in terms of the nodal dis- 
placements 

{a} = [A](u] (114) 

where 





This leads to an expression of strain in terms of the independent 
degrees of freedom (nodal displacements), 

[e] = [WHu] (115) 

where [W] = [B][A], 
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4.6.3 


Stiffness Matrix 


Integrating the elastic strain energy over the whole ring 
e lemen t 


U = i 


{a}Me)dv 


V 


(116) 


where {e} can be expressed as the difference between the total 
strain and initial strain 

{€} = - { 6 ^} 


and using the principle of virtual work we arrive at the stiff- 
ness matrix 


[k] = 2tt 


[W]'^[C] [Wjrdrdz 


(117) 


V 


where [C] is the matrix of elastic constants (see Section 


4.6.7) . 

In order to integrate Eq. (117) we must integrate the fol- 
lowing quantities 


— drdz = I, 
r 1 


^ drdz = I^ (118) 


2 

“ drdz = I^ 
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Carrying out this integration we arrive at 


- log 


^i + 


(A 


kj 


- A. .) 
Ji 


log r. + 


^^ik 


A, .) 


log r. 


+ B..(r. - r.) +B.,(r. - r, ) + B.,(r, - r.) 

ij X j' jk' j k"^ ik^ k 


lo = 


((A?^- A^j^) log r^+ (A^j- Aj^) log r^ + (A^j^ - A^^ ) log rj^J/2 

+ A . 1 B . 1 (r, - r . ) + A, . B, . (r . - r, ) + A . . B . . (r . - r . ) 

ik xk' k x' kjkj'j k' jijx'i j' 


, r^2 fl , Jl f 1 2. , ^2 . 2 2.,,, 

+ f®ik<’^k - + Bkj<’^j - >^k> + 


I. = 




2 2 2 

■*■ ^ik^ik^’^k ' ’^i^ ■*■ " ’^k^ ^ji^ji^’^i " 

+ 'Aik®L(4 - 4 ^ + \j®kj<’^j - 4> + 


+ tBj^(r3 . ^ 33_(^3 . ^3^ ^ ^3^(^3 . ^3^,/, 


Two special cases arise. One when r^ = 0 which affects the 
terms involving (Aj^ ~ A^j) log r^. Using L'Hopital's rule the 
limit of these terms is zero and can be taken into account by 
omitting the logarithmic terms from I^, I 2 , and . The 
second special case arises when r^ = r^^ in which case we do 
not integrate over the area under i-k. This can be accomplished 
by setting A^j^ = B^j^ =0 in I 2 , and These points 

are covered in Ref. 45. Note that in the expression for 


in vjhich case 


except for the case in which r^. = 

ij'-l - 3 ' ^ "ik^"k 


1 

B,.(r, - r.) + B,, (r, - r,) + B,, (r, - r, ) = z, - z. 


4.6.4 Initial Strain Matrix 

The plastic load vector is given by 


fQ] = 2 tt 


[W]^ rdrdz [C]{€^] = [k*]{e^] 


(119) 


V 


where [k ] is the initial strain stiffness matrix and is 

the vector of plastic strains taken to be constant within the 
element and evaluated at the centroid. 

The integration associated with the initial strain matrix, 
Eq. (119), can be found by evaluating the integrand at the cen- 
troid, i.e.. 


[Q] = 2v 


[W] rdrdz [C][e^} 


{Q] = 27rFA[W]'^[C][e^} 

where A is the cross sectional area, F is the centroidal 
point, i.e., F = (r. + r. + r,)/3 and [W]"^ is the [W]"^ 

1 J K 

matrix evaluated at the centroid. 


( 120 ) 


4.6.5 Stress Calculation and Evaluation of Plastic Strains 

Stresses and plastic strains are calculated at the centroid 
of the element. All calculated quantities such as stress, plastic 
strain, shift in yield surface, etc., are output in the principal 
direction of orthotropy as defined by the angle p (see Fig. 19). 
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4.6.6 


Thermal Stress Calculations 


Orthotropic thermal stress calculations are allovjed. Dif- 
ferent thermal coefficients of expansion may be input in the 
principal directions of orthotropy. Temperatures are input at 
nodes. Since a constant thermal strain distribution is assumed 
within the element, the temperatures at each of the three nodes 
are averaged. This average value is assumed to hold throughout 
the element. The thermal load vector is obtained by multiplying 
the thermal strains by the initial strain stiffness matrix 


(APt,) = [k-HA.^,) 


where 




= 




QsT 


0 

0 

0 
V 


4.6.7 Material Properties 


( 121 ) 


We assume an orthotropic material whose principal directions 
(1, 2,3) correspond to (1, 2, r) as shown in Fig. 19. We can 
then express the stress-strain relationship in the principal di- 
rections 

Ca'} = [C][€'} (122) 

where corresponds to the principal direction whose components 
are 111, 22, rr, r2, rl, 12 j and 
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^11 

= Ej^(l 

''32^23 

^)/A 


^12 

= C21 = 

Ei(Vi2 

+ ^^ 32 ' 

'13)/^ 

Ci 3 

= Sl = 

Ei(Vi3 

+ v^ 2 ^ 

'23) M 

^22 

= E2<1 

''l 3''31 

)/A 


^23 

" ^32 

® 2^''23 

+ ^13' 

-21) M 

C33 

= E3(1 




•^44 " ^23 
^55 " ^13 
^66 " ^12 




A = 1 - 

^'^12^23^31 “ ^ 

13^31 

^12 

^21 

V V 

23 32 


We no\^ 

express the 

stresses and 

strains in 

the global coordinates 

as 



[a] = 

[CHe 

} 



(123) 

where 









[C] = 

[R] 

’^[C][R] 









2 

cos p 

sin p 

0 

0 

0 

-sinp 

cosp 



2 

sin p 

2 

cos P 

0 

0 

0 

sinp 

cosp 



0 

0 

1 

0 

0 

0 


[R] = 


0 

0 

0 

cosP 

sinp 

0 




0 

0 

0 

-sinp 

cosp 

0 



.2 

sinP cosp 

-2 sinp cosp 

0 

0 

0 

2 

(cos P 

- sin^P) , 
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In the current version of the REVBY program for strain 
hardening materials only initially isotropic material behavior 
can be specified if a plastic analysis is desired. Either a 
linear or nonlinear hardening law may be chosen. For ideally 
plastic materials, isotropic or orthotropic plastic material 
properties may be specified. At the end of each half -cycle the 
plastic material properties may be changed. 

The thermal coefficients of expansion may be specified to 
be orthotropic in the principal directions. These may not be 
changed at the end of each half -load cycle. 


4.6.8 


Loads 


Concentrated and Line Loads — External nodal forces 
are expressed in force per unit length of circumference except 
at 

r = 0 there must be F = F„ = 0, and 

IT C7 

force, in pounds for example. 


r = 0. The total load applied is therefore 2Trr{F^} . At 


F^ represents actual 


Surface Tractions — Distributed surface force/unit area can 
be specified at nodes in the r, z, and 6 directions. If we 
assume a linear variation of traction across the face, i.e., 
p = a + bi (Fig. 21), we arrive at the consistent load vector 


(124) 


— > 

F. 


1 1' 


"p-1 

X 


3 6 




} = Tr(r^+r^)L 

1 1 

< 

1 

! 

1 1 

F. 


.6 3 




where 


p^ corresponds to the traction at node 


1 . 


Note: All loads are incremented proportionally from the critical 

load . 
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Fig. 21 Distributed Surface Forces 


4.6.9 General Comments 

• Note that for orthotropic material properties, the R-Z 
axes and the axes of principal orthotropy must be coincident, 
i.e., no rotation about the 6-axis is allowed. 

• This element is available in the REVBY module. 

4 . 7 Membrane Triangles 
4.7.1 Introduction 

The membrane triangles described in this section are classi- 
fied into the following three categories: 1) 3-node constant 

strain, 2) 6-node linear strain, and 3) 4- and 5-node hybrid 

elements. The terms constant, linear, and hybrid refer to the 
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strain distributions that exist in the element as a consequence 
of choosing an assumed displacement variation. A brief descrip- 
tion of the elements associated with these three categories 
follows . 

Constant Strain Triangle (CST) — This well-known plane 
stress membrane element is one of the most widely used elements 
for the idealization of membrane structures. Its derivation is 
based on the assumption of a linear distribution for the in-plane 
displacements u and v, and consequently, leads to a constant 
strain state within the element. Each vertex is allowed two de- 
grees of freedom (the in-plane displacements u and v) for a 
total of six degrees of freedom for the element. Consistent with 
the total strain distribution, the initial strains (plastic 
strains) are assumed to be constant within each element. Stiff- 
ness and initial strain matrices have been developed and success- 
fully used in Refs. 6, 8, and 30. 

Linear Strain Triangle (LST) (Fig. 22) - In regions of high 
strain gradient, the CST triangle is not sufficiently accurate to 
be used in a plasticity analysis unless a very fine grid is em- 
ployed. The linear strain triangle (LST) remedies this short- 
coming of the CST element. The assumption of a quadratic dis- 
tribution for the in-plane displacements allows for a linear 
strain variation within the triangle. Two degrees of freedom at 
each node (u,v) for each of the six nodes (three vertex and 
three midside nodes) give this element a total of 12 degrees of 
freedom. The initial strains are assumed to be constant within 
each element and are evaluated at the centroid. Both stiffness 
and initial strain matrices have been developed and successfully 
used in Refs. 8 and 30. 
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Fig. 22 Linear Strain Triangle 

Hybrid Triangles (HST) — In transition regions, i.e., re- 
gions in which stresses and strains change from rapidly varying 
to slowly varying, it becomes convenient and efficient to switch 
from linear strain triangles to constant strain triangles. This 
is accomplished by using four-and five -node triangles to main- 
tain compatibility with both the CST and LST elements. These 
elements together with the CST and LST elements were originally 
used in Ref. 4 and are referred to as the TRIM 3 through TRIM 6 
family. For these mixed formulation hybrid elements, the dis- 
placements along edges may vary quadratically or linearly, de- 
pending on whether an LST or CST triangle is contiguous to the 
respective sides. Again, the plastic strain distribution is 
assumed constant within each element. 
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4.7.2 


Displacement Assumptions 


Using a local centroidal coordinate system with the x-axis 
parallel to side 1-2 shown in Fig. 22, the in-plane displace- 
ments u and V may be written in terms of the area coordinates 
0^1 , 0^2 > • 

U = CDj^( 203^ - l)u^ + g:)2(2(X>2 " “ ^^^3 

(125) 

+ 4co2^i>2'^4 ^“3*^1^ 5 ^^1^2 6 

or in matrix form as u = [N]{u^} where [N] is a matrix of 
shape functions representing the assumed displacement field 
within the element. 


The terms u^ (i = 1-6) refer to nodal displacements at the 
vertex and midside nodes, the subscript o refers to the nodal 
displacements, and the area coordinates are written in tertns of 
the local x,y coordinates in the following matrix form 

I'' 




1 ^ 2 )^ A 




^2^3 ’^3^2 


X3yi-Xiy3 

Xiy2-X2yi 


72-73 


73-71 


7i-72 


X 3 -X 2 


X 1 -X 3 


X 2 -X 1 
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(126) 


or 


ri> 


2A^ 

*>1 


CM 

3 

1 

> 2A 

2 A 2 

^2 

^2 

l“3. 


2A^ 

*>3 

^3. 


X 


y) 


An expression similar to Eq. (125) can be written for the v 
component of displacement. The particular form of Eq. (125) 
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represents a quadratic poljmomial in tlie coordinate variables x 
and y and consequently is associated with the six -node LST 
element. A linear pol 5 momial representation of the displace- 
ments (CST element) is obtained from Eq. (125) if we enforce the 
following conditions 

The hybrid (HST) element is obtained by enforcing any one, or 
two, of the above conditions. 


4.7.3 Stiffness Matrix 


Neglecting the effects of large geometry changes, the strain 
displacement relations are determined by applying the appropriate 
differential operator on the displacement functions as 


X 


\A. 

dx 

0 

/ e 

' y 

II 

II 

0 

dy 

T 






dx. 


The strain-nodal displacement relations may 
in matrix form as 


(127) 


be represented 


(e} - [D][NHu^} = [W][u^] 


(128) 


For the LST and HST element the [W] array is a linear com- 
bination of the coordinate variables x and y and may be 
written as the sum of the following terms 


[W] = [W^] + [W2]x + [W3ly 
3em 3«m 3*m 3*m 


(129) 
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where m may be 8, 10, or 12. Note that the [W] array will 
vary with the number of nodes used for the element. For example, 
for the four-node element, appropriate transformations associ- 
ated with the displacements removed from two of the three midside 
nodes must be applied to the full (3x12) [W] array to obtain 

the correspondingly reduced (3x8) array. Correspondingly, 
the strain-nodal displacement equation for the CST element is 
represented by an expression of the form of Eq. (127) with 

[W] = [W^] (130) 

3«6 

The local stiffness matrix for an element of uniform thick- 
ness is represented in the following form 


[k] = h 


''A' 


[W]'^[E ] [W]dA 


(131) 


where [E] is the array of elastic coefficients relating stress 
and strain components in local coordinates. 

For the LST and HST elements the form for [W] given in 
Eq. (129) leads to the following expression for the local stiff- 
ness matrix 


[k] = Ah[W^]^[E][W^] + yI (xj + X2 + X3)[W2]^[E][W2] 
m#m 

+ ^2^2 ^3y3^ '■^3^ [W3]^[E] [W 2 ]) 

■*■ 1 ! ^^1 ^2 y ^) [W3]"^[E] [W3] 


(132) 


where m = 8, 10. or 12. Similarly, for the CST element the 
expression for [k] reduces to 
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(133) 


[k] = Ah[W ]’'’[£] LW ] 

6.6 ° ° 

The local stiffness matrices [k] are then transformed to 

T 

global coordinates by premultiplying them by [T] and post- 
multiplying by [T], where [T] transforms the displacements 
from local to global quantities. 


4.7.4 Geometric Stiffness Matrix 


For two dimensional planar problems it is conjectured that 
large geometry changes (e.g., localized near a crack tip) are 
basically those involving large rotations. Additionally, the 
presence of large membrane stresses significantly affects the 
out-of-plane stiffness of these elements when they are used for 
built-up structures. This is important in analyzing large de- 
flection and stability problems. On this basis, the inclusion 
of small strain, large rotation effects (geometric nonlinearity) 
is effected through the use of the following nonlinear incre- 
mental strain-displacement equations 


Ae 

X 

Ae 

y 



(134) 


A _ 5 , 4u oAv ftAw P^Aw 

'^xy ~ dy dx dx dy 


where A represents an incremental change between two adjacent 
states of the loading history, and 


CD 

z 


- l( ^ ^ 

^\dx dy/ 


(135) 


130 


The displacement w is in the local z-direction. The dis- 
placement assumption used for w is the same as that used for u 
and V in corresponding elements. 

The nonlinear terms of Eq. (134) give rise to the geometric 
(or initial stress) stiffness matrix. Its local form may be 
represented as follows for an element of uniform thickness 




[n] ^[2][n]dA 


(136) 




The matrix [fi] is used to express the rotations in terms of 
nodal displacements, e.g.. 


' Aw 

/ Aw 



[fiHAu^} 




(137) 


where includes the w terms and [ 2 ] is, using a 

convected coordinate (updated geometry) approach for the solu- 

tionof geometric nonlinear problems merely an array of stress 

resultants N , N , N , for each element (or node) . 

X y xy 

Centroidal values of the stress resultants are used in all 
order elements. The local geometric stiffness matrices of 
Eq. (136) are then transformed into the global system by pre- 
multiplying them by [T] and postmultiplying them by [T]. 


4.7.5 Initial Strain Stiffness Matrix 

The treatment of plasticity within the PLANS system requires 
the incorporation of two basic assumptions. The first is, in 
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determining the effects of plastic flow within a body, plastic 
strains are interpreted as initial strains; the second is that 
total strains (or their increments) may be decomposed into 
elastic and plastic components. The local initial strain matrix 
that appears in the equilibrium equations may be written in 
terms of previous variables as 


[k*] 


= h 


[W]'^[E] [N*]dA 


A 


(138) 


where [N ] represents the matrix of shape functions used to de- 
scribe the assumed plastic strain variation within an element. 

If we choose to consider the plastic strain state at the centroid 
of the element as representative for the element (regardless of 
the total strain variation) then Eq. (138) may be written as 

[k*] = Ah[W, ][E] (139) 

m®3 


where m = 8, 10, or 12 for the LST and HST elements, and 

[k*] = Ah[W^][E] (140) 

6e3 ° 


for the GST element. 

Plastic Load Vector — The product of the initial strain 
matrix and the vector of plastic strains represents the "effec- 
tive” plastic load for the structure. This load is applied in 
addition to the prescribed set of mechanical loads to determine 
the displacemeiit state. The aforementioned product is determined 
at the element level; that is, rather than assembling an initial 
strain coefficient matrix for the entire structure and multiply- 
ing by the vector of plastic strains, the vector product of the 
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element initial strain matrix and plastic strains for the element 
is assembled for the entire structure. Thus, the local plastic 
load vector may be represented in the following form 


X 


N 


(Q) - I 

i=l m*3 


(141) 




where m = 6, 8, 10, or 12, N = number of plastic elements, 


and 


eP, tP, 

x’ y 


and 7 ^ represent plastic strain components. 


xy 


This vector is then premultiplied by [T] 
global coordinates. 


to transform it to 


4.7.6 Stress Calculations 

The total strains for each element are evaluated at the 
centroid using the strain-displacement relations previously dis- 
cussed. For small deformations the linear relation becomes 

[e] = [W^Ku^} (142) 

3«m 

where m = 8, 10, or 12 for the LST and HST elements, and 

(e} = [W^][u^] (143) 

3*6 

for the C'ST element. Once again, note that [W^] will vary de- 
pending on the value of m. 

The nonlinear contribution to the strains, as expressed in 
Eq. (134) requires the array to be evaluated at the cen- 

troid. Stresses and plastic strains are also evaluated at the 
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centroid for LST, HST, and CST elements. Although the stresses 
in an LST or HST element are not uniform, but vary from node to 
node within an element, the centroidal value is used to repre- 
sent the stress state of the element. The plastic strains, on 
the other hand, are a priori assumed to be constant within the 
element regardless of the element type (e.g., LST, HST, or CST). 
It should be noted that using a constant centroidal value of 
stress for the LST and HST elements does not degrade the stiff- 
ness influence coefficients for these elements.' It does, how- 
ever, influence the nonlinear response of the structure in that 
plastic strains are determined from these stresses. The inaccu- 
racies associated with using centroidal stresses is most pro- 
nounced in regions of rapid stress gradients. 

4.7.7 Thermal Stress Calculations 

Nodal temperature input is converted to average element tem- 
peratures and subsequently to centroidal thermal strains for LST, 
HST, and CST elements. Orthotropic thermal coefficients of ex- 
pansion may be specified in the principal material directions of 
orthotropy. The thermal load vector is determined in the same 
manner as the plastic load vector, e.g., as the product of the 
initial strain matrices and the vector of thermal strains for the 
element and summed over all the elements of the structure. 

4.7.8 Material Properties 

All elastic material properties are constant within the ele- 
ment and are assumed to be temperature independent. Orthotropic 
properties may be specified by defining the necessary material 
constants and the direction of the principal axis of orthotropy, 
P, for each element. In general the stress-strain relation may 
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be written in an arbitrarily rotated, orthogonal, 1-2 plane 
(Fig. 22) as 


where 


<’11' 

1 

^11 

*^12 

0 ■ 

rill 

'i 

^22^ 

^ = 

^21 

^22 

0 

( ^22 ) 

Jl2. 
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^33. 



11 

= Ei/(1 - 


12 

= "12^11 
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= E^/Cl - 

'' 21 ' 12^ 

21 

" ^21^22 


33 

= Gj_2 



(144) 


2 

and ^11^22 ^12 ^ ^1’ ^2’ ^12 ^ ^ for a positive defi- 

nite stiffness. 


Thermal coefficients of expansion in the 1 and 2 direc- 
tions may be specified independently. 

The von Mises yield criterion is used to define the limits 
of elastic behavior for isotropic materials. Hill's yield condi- 
tion is used for materials having different yield stresses in the 
two orthogonal global directions. This yield criterion for plane 
stress is written as 

f = (G+H)o^ + (F+H)cJ 22 ~ 2Hct^^22 2Nt^2 ” ^ 


where 
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G+H = 1/X^ 

F+H = 1/Y^ 

2H = 1/y? + 1/Y^ - 1/Z^ 

2N = 1/T^ 

and X,Y,Z represent the normal tensile yield stresses in the 
1, 2, and 3 directions, respectively, and T is the yield 
stress in shear in the 1-2 plane. Hill's condition, Eq. (134), 
reduces to the von Mises yield condition for isotropic yield 
stresses. Furthermore, the following condition must be enforced 
to ensure the convexity of the orthotropic yield surface with 
respect to the origin in stress space 

2 

^ + 4 ' ’t) >0 (146) 

X Y X Y Z 

Three options are available to describe the nonlinear mate- 
rial behavior: 1) elastic-ideally plastic, 2) elastic-linear 

hardening, and 3) elastic-nonlinear hardening. For the first 
option only the yield stresses, X,Y,Z, and T need be speci- 
fied. For hardening behavior the kinematic hardening theory as 
proposed by Prager (Ref. 42) and modified by Ziegler (Ref. 37) is 
used to describe the subsequent hardening behavior of the mate- 
rial beyond initial yielding. For linear hardening the slope of 
the tensile stress-plastic strain curve is specified (tangent 
modulus) . If nonlinear hardening is desired the Ramberg- Osgood 
parameters (Ref. 38) must be specified (see Section 3.4 for fur- 
ther discussion). For orthotropic kinematic hardening the 
Ramberg-Osgood parameters must be specified for each of the 
three stress components. 
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4.7.9 


Loads 


The following mechanical loads may be applied to the mem- 
brane elements: 

• Concentrated forces and moments applied at specified 
nodes in the global directions in the units of force 
or force times length. 

• Distributed edge loads in the plane of the element. 
These are assumed to vary linearly from node to 
adjacent node along the edge on which the load is 
applied. The magnitudes of the loads are specified 
in the global directions at each pair of adjacent 
nodes . 

The consistent load vector, [P], associated with the dis- 
tributed load case is determined from the shape function, [N], 
Eq. (125), 


{P] 


[Nj'^’CpldA 




(147) 


4.7.10 General Comments 

• The same assumed constant plastic strain distribution is 
used for all of the elements of the membrane type (CST, HST, and 
LST) . Although this distribution is compatible with the total 
strains in the CST elements only, it should not introduce a sig- 
nificant degradation of the accuracy associated with the HST and 
LST elements. 

• Centroidal values of stresses and strains are used for all 
the elements. 
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• Although the situation illustrated in Fig. 23 is allowed 
(no error message triggered) it is not recommended. The dis- 
placement field along the common edge A-B from each of the two 
adjacent elements will not be compatible. 

• The applicable modules for these elements are BEND, OUT-OF- 
PLANE, and OUT -OF -PLANE -MG. They are also available in a special 
fracture module called FAST as well as in preliminary versions of 
OUT-OF-PLANE designated as PLANE. 


quadratic edge displacement 



Fig. 23 Incompatible Edge Displacements Along A-B 
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4.8.1 


Introduction 


This element is used for analyzing thin-v^alled structures 
where local bending as well as membrane effects are important. 
The element is of constant thickness, triangular planform, and 
has 12 degrees of freedom at each node (see Fig. 24). It is 
obtained by superposing the Quadratic Strain Membrane Triangle 
of Felippa (Ref. 16) and Tocher and Hartz (Ref. 48), with the 


y , V 



Fig. 24 Triangular Plate Element - Bending and Membrane 



quintic displacement bending triangle of Bell (Ref. 47) and 
Cowper et al. (Ref. 50). The local coordinate system is defined 
at the centroid of the element. The local x-axis is parallel 
to the 1-2 edge. The local y-axis is perpendicular to the 
x-axis and in the plane of the element. The local z-axis is 
perpendicular to both of these axes. 

4.8.2 Displacement Assumptions 

, , , , 3 

u = a^ + a^x + a2y + . . . + ^-^qY 

V = b^ + b^x + b2y + ... + b^Qy^ (148) 

3 5 

vj = c^ + c^x + c^y + ... + c^^y + ... C2^y 


The in-plane displacements u,v are complete cubic poly- 


nomials. The 20 undetermined coefficients (a^ and 
related to nodal degrees of freedom through the matrix 
These nodal quantities are u. , v. , u , v , u 

1. 1 « • A "'V • 

’ 1 ’ 1 ’*^1 

i = 1, 2, 3, and centroidal displacements 


bi) 


are 

-1 


V 


where 


= CAJ 


-1 


m 


(149) 


The 20 in-plane degrees of freedom associated with u and 

V are reduced to 18 in-plane vertex nodal degrees of freedom 

by Gaussian elimination (equivalent to static condensation) of 

the centroidal degrees of freedom u^,v^. The remaining nodal 

degrees of freedom are u. , v. , 0 , € , € , e , where 0 
° r’l’z.x.y.xy. z. 

1 X -^1 -^1 1 

is the local rotation about the z-axis, 0 = ^(-c'u./ciy + 

1 

Sv./c*x), and e = Su./c)x, e = ^v./c)y, e = 4(c^u./t^y + 

JL <4^ • -L y • 0 JL 

X X X 

^v^/5x) are the tensor components of the membrane strains. 
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The out-of-plane displacement w is a complete quintic 
pol 3 nnomial in x,y. The 21 coefficients (c^) are related to 
the vertex nodal degrees of freedom w., v; , v; 


> w 

jXf ,y^’ ,xx^ 


w 


w 


.yXi >xy^ 


, and the normal midside nodal rotations vj 


through the matrix 


-1 




- 1 , 


ta^) = (A^l ^[d^) 


(150) 


The 21 degrees of freedom are reduced to 18 by requiring 
that the normal slopes vary cubically along an edge (see Ref. 49), 
thus eliminating the three normal slopes at the midside nodes. 

The local bending degrees of freedom remaining at each vertex 


node are , 6 = w ,9 = -w , and k = w , k 

1 X. ,y, y. ,x, X. ,xx. y, 

1 ■'i ■' 1 . 1 X 1 -^i 


w 


K 




= V? 




The bending and membrane stiffnesses are then superposed to 

obtain a 36 degrees of freedom element vjith the following 

12 local degrees of freedom at each node u, v, w, 9 , 6 , 0 , 

O 5 J 5 y 5 2 ’ 

e,€,6 ,K,K,K 

X y xy X y xy 


4.8.3 Formation of Stiffness Matrix 

The membrane and bending stiffness matrices may be written 


as 


[k ] = h[A ] 
m m 


[W ]’^[E][W^]dA[A^]'^ 
tn m m 


[kbl [C]^[A,]-T 


[W^]'^[E][W^]dA[A^]"^[C] 


(151) 
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Here h is the plate thickness; [W ] and [W, ] are the 

mb 

matrices relating the strains to nodal parameters as 

^ ^ ] 

X 1 1 dx 


m ' 


= / 




5u 

c)u j 

5y r^x ' 


= [W^Ha } 

m m 


[k^] = ^ K 


( w > 

, XX ( 


' ( «,yy ■'= 


(152) 


\ 


2v/ 


,xy; 


and [C] is the condensation matrix for the bending component. 
[E], the matrix of elastic constants relating stress to strain, 
is discussed more fully in the section on material properties. 
After formation of the full 20 degrees of freedom membrane 
stiffness matrix it is reduced and reordered. The area integra- 
tion for bending and membrane is performed using Gauss-Legendre 
integration of fourth order. 


4.8.4 Formation of Geometric Stiffness Matrix 

The geometric stiffness matrix was developed to be used with 
an updated coordinate system approach (Ref. 26) for solving geo- 
metric nonlinear problems. Small-strain moderate rotation as- 
sumptions were made to arrive at the strain-displacement equa- 
tions used. These are 
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where cjd , cn , and co are Infinitesimal rotations, 

X y ’ z 

The component of the geometric stiffness that affects 
bending stiffness is written as 


[kcb] = [C]’^[Ab]-^ 


[w, ]'^[N][W, ] dA[A, ]"^[C] 


Here 


[N] = 


N 


and 


( w 
1 ,x 
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. N 

I xy 
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xy 


N 




where relates the slopes to the nodal parameters a^ 

geometric stiffness that affects the in-plane stiffness is 

-T I rr', 1"! 


[k_ ] = (N +N ) [A ] 
Gm ^ X y' m 


[W„]^[W^]dA[A^]' 
mm m 


(153) 


the 


( 154 ) 


. The 


( 155 ) 
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A A 

where co = [W ]{a } and [W ] is a row vector relating the 
z m m m ^ 

rotation to the nodal parameters 

These matrices are added to the conventional stiffness 
matrices and reordered and reduced in the manner discussed in the 
description of the conventional stiffness matrices. For problems 
involving geometric nonlinearity all stiffness matrix integrations 
are performed using fifth order Gauss-Legendre integration over 
the triangular regions. The values of N^, N^, and used 

to form the geometric stiffnesses are obtained by averaging re- 
spective nodal quantities from each of the three nodes. 


4.8.5 Formation of Initial Strain Matrix (Plastic Load Vector) 

The membrane and bending initial strain matrices are written 
as (Ref. 26) 


m 


m 


-T 


[W^]^[E][C]dA 




[k*] = [C]^[A^]"'^ 


[W^]^[E][C]dA 


A 


(156) 


The initial plastic strains are assumed to vary linearly 
from node to node within the element and have an arbitrary varia- 
tion through the thickness at the nodes, i.e., 

{e^(x,y,z)l = [C(x,y) ]i;e?(z)} (157) 

where [C(x,y)] is a linear function matrix in (x,y) relating 
the plastic strains in the element to the nodal plastic strains. 
The area integration for the bending and membrane components of 


144 



r 


the initial strain stiffness matrices is performed using Gauss- 
Legendre integration of fourth order. The membrane component of 
the initial strain matrix is condensed from 20 to 18 degrees 
of freedom using Gaussian elimination. 

The plastic load vector CAQ} is the product of the initial 
strain matrix and the vector of initial strains or, for plastic 
analysis of plates (and shells), inelastic forces and moments, 

{AQ} = [k*HAe} (158) 


where 


(Ae} = 




/ ^ dz 
f ^ dz 

J y 

r Av dz 
•' xy 

f Ae zdz 
f Ae zdz 

J y 

f AyP zdz 

J jry 


— 


/Ae dz 

•' X 


rAyf’ zdz 

j xy 




f AyPyZdz 


The plastic strains are evaluated at the nodes of the tri- 
angle and integrated through the thickness using Simpson's rule. 
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Up to 20 layers (21 points) are alloT^ed. The number of layers 
must be an even number. If the number of layers is not speci- 
fied then the default of 10 layers is used. This has been 
found to be sufficient for most problems. 

4.8.6 Stress Calculations and Evaluation of Plastic Strains 

Stresses are calculated at each of the three nodes of the 
triangle. There is no averaging of membrane strains or curva- 
tures from adjacent triangles. Since the orthotropic material 
properties are given in principal directions of anisotropy, all 
calculated stresses are in the element material coordinate sys- 
tem defined by the angle P (see Fig. 24). These directions 
must also be used in conjunction with Hill's orthotropic yield 
criterion. Stresses are also calculated at each integration 
point through the thickness and will be printed out at these 
points if requested. Otherwise just top and bottom surface 
stresses and strains are output. 

4.8.7 Thermal Stress Calculations 

Orthotropic thermal stress calculations are allowed. Dif- 
ferent thermal coefficients of expansion in the principal direc- 
tions of orthotropy may be specified. A parabolic temperature 
distribution through the thickness is assumed. Temperatures at 
the top, bottom, and middle surface through the thickness at 
each node are input. Like the plastic strains, temperatures are 
assumed to vary linearly from node to node. Temperatures at each 
layer through the thickness are determined from the assumed para- 
bolic distribution. The thermal load vector is obtained by 
multiplying the thermal strains by the initial strain stiffness 
matrix 
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The integrations are evaluated in closed form based on the para- 
bolic temperature distribution. 


4,8.8 Material Properties 

All material properties are constant within the element and 
are assumed temperature independent. Elastic properties are as- 
sumed to be orthotropic with the direction of the principal axis 
of orthotropy, P, specified for each element 
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2 

where ^n^22”^12 ^ ^l’^2’*^12 ^ ^ positive definite 

stiffness. Thermal coefficients of expansion in the 1 and 2 
directions may also be specified independently. 
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Plastic material properties may be assumed to be ortho- 
tropic with Hill's yield criterion for plane stress used for 
initial yielding. This yield criterion is written as 

f = (G+H)aJ + (F+H)02 - 2Eo-^02 ^ 

where 

G+H = \ 

X 

F+H = ^ 

Y 



and X,Y, and Z represent the normal tensile yield stresses in 
the 1,2, and 3 directions, respectively, and T is the yield 
stress in shear in the 1-2 plane. The following additional 
stability criterion must be satisfied by the yield stresses for 
Hill's equations to be used 





(162) 


Kinematic hardening theory as proposed by Prager (Ref. 42) and 
modified by Ziegler (Ref, 37) is used to describe the subsequent 
hardening behavior of the material. 

Three different plasticity options are available: 1) elas- 

tic, ideally plastic, 2) elastic linear strain hardening, and 
3) elastic nonlinear strain hardening using a Ramberg-Osgood 
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power law representation of the stress-strain law. For ideally 
plastic materials only the yield stresses X,Y,Z,T need to be 
specified. For linear strain hardening the slopes of the 
plastic portion of the stress-strain curves (tangent moduli) 
must be input for each of the three stress components. If non- 
linear hardening is desired then two Ramberg-Osgood shape parame- 
ters (see Section 3 for further discussion) for each stress com- 
ponent must be specified. 

Different stress-strain components may have different 
hardening laws, e.g., hardening, 

linear hardening, and '’^xy’^xy hardening. However, 

if one component is elastic-ideally plastic all must be. Addi- 
tionally, at the end of each half -cycle of loading, the plastic 
material properties may be changed. Elastic properties may not 
be changed at the end of each half-cycle. 

4.8.9 Loadings 

The following mechanical loads may be applied to the plate 
element: 

« Concentrated forces and moments 

® Distributed edge loads perpendicular to and 
in the element plane 

® Distributed surface loads perpendicular to 

and in the plane of the element. 

Concentrated Forces and Moments — These are applied at 
specified nodes in the global directions. 

Edge Loads — These are assumed to "vary linearly from node to 
adjacent node along the edge on which they are applied. They are 
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applied in the element local coordinate system. The three com- 

th 

ponents of the distributed edge load are specified at the i 
th 

and j node of the edge. Edge 1 is defined to be the side of 
the element connecting nodes 1 and 2; edge 2 is defined to be 
the side connecting nodes 2 and 3; and edge 3 is defined to be 
the side connecting nodes 3 and 1 (Fig, 24). Applicable members 
with the same distribution along the respective edges are then 
specified. Consistent load vectors for these loadings are cal- 
culated in the element routines 
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These are then reordered and reduced as described in the initial 
strain section. 


Distributed Surface Loads — The three components of the sur- 
face loads are assumed to vary linearly from node to node in the 
plane of the element. They are applied in the element local co- 
ordinate system perpendicular to and in the x,y directions. 

The three components of the surface load P^jP^jP^, ^tre specified 
at the lj2,3 nodes of the element. Applicable members with the 
same values of components of the respective nodes are then speci- 
fied. Consistent load vectors are calculated in the element 
routines. These are determined from the following relations. 
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The load vectors can be -written as 
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These are then reordered and reduced as in the initial strain 
section. 


4.8.10 Equilibrium Correction 

We may write the equilibrium correction terms for the plate 
element as (Ref. 36) 


where 
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Assuming a linear stress variation from node to node this can be 
rewritten as 
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The component is formed in the load vector formation 

and merely consists of multiplying the consistent load vector by 
the multiplicative scalar appropriate to the current load value. 
The second component is virtually identical to the initial strain 
matrix and is formed at the same time. It is multiplied by the 
vector of stress and moment resultants. The latter is obtained 
by numerically integrating the stresses through the thickness 
using Simpson's rule. 


4.8.11 General Comments 

• In general, displacements and rotations between elements 
along an edge will not be compatible if the adjacent elements 
are not in the same plane. This is due to the assumption of a 
quintic polynomial representation for the out-of -plane displace- 
ment w and a cubic representation for u and v. If the 
structure is planar then displacements and rotations will be com- 
patible and membrane strains and curvatures will be identical in 
elements with common nodes. 

® The displacements and rotations are assembled in the global 
system. Membrane strain and curvature degrees of freedom are as- 
sembled in a "local global system." The local coordinate system 
of the first element specified that contains a node becomes the 
reference system for that node. For this reason the out-of-plane 
angles between elements should not be large and each grouping of 
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elements about a node should be a "shallov? shell." If this is 
not so, for example, at a node on the corner of a box, the mem- 
brane strains and curvatures should not be assembled with those 
of an element with which there is a large angle. To avoid link 
ing all of the degrees of freedom of the two nodes, multiple 
definition of a point with two node numbers can be made and the 
same coordinates assigned to the two nodes. Multipoint con- 
straints can then be used to selectively link displacement and 
rotation degrees of freedom to ensure continuity of these quan- 
tities . 

• This element is available in the BEND module. 


4 . 9 Isoparametric Hexahedron 
4.9.1 Introduction 


A review of the present state of isoparametric elements is 
presented in Ref. 51. These elements are characterized by a 
mapping procedure used to map a simple shape in the local or 
natural system into a curvilinear shape (actual shape) in the 
global system. For the three dimensional solid- element pre- 
sented here, a cube in local coordinates is mapped into a hexa- 
hedron in global coordinates (Fig. 25). The mapping function can 
be written in the form 


\y r ^i) 


(167) 


where x, y, and z represent the global coordinates, ^,C, 
and T) represent the local coordinates, and x^ represent the 
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a) Local Coordinates 

Fig. 25 


b) Global Coordinates 


Isoparametric Hexahedron 


nodal coordinates in the global system. The mapping function, 
or shape functions,' N^, have the property of being unity at 
node i and zero at all other nodes. 

4.9.2 Displacement Assumptions 

The element deformation is described using the same shape 
functions as in Ref. 51, hence the name isoparametric 

^ N^(^,C,Ti)<t . (168) 

i 

where 4)^ represents the unknown generalized degrees of freedom 
associated with an individual element. 

With the use of "relative coordinates" we now derive a 
single isoparametric hexahedron which has a variable number of 
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nodes, of between 8 and 20. In practice this is especially 
useful where it is necessary to generate a mesh containing adja- 
cent elements with different numbers of nodes, as occurs when a 
mesh changes character, e.g., going from a more coarse to a 
finer mesh. For a more detailed description of the use of rela- 
tive coordinates see Ref. 52. 

For an isoparametric hexahedron of between 8 and 20 nodes 
the shape function is as follows 

<t> = ^ ^i'*’i ^ •••> S plus any midside nodes (169) 

i 

where for the midside nodes (i = 9, ..., 20), as they 

exist, and nT = P. - ^(Pt + P, + P,,) for the corner nodes 

1 X I J K 

(i = 1, . . . , 8) where I, J, and K represent the midside 

nodes, as they exist, adjacent to corner node i. 

The P functions are given as follows 

Pi = id + e„)(i + C„)(1 + -n^) (170) 

for the corner nodes and 

Pi = id - e^)d + C^HI + ■ o''i 

Pi = i(i + - ?^)(1 + Oo) <^1 - 0 ;> d7i) 

Pi - i(i + e^Hi + - 0^) ’ll - 0^ 

for the midside nodes, where ^ H and ri = titi . . 

’ o l o X 'o *'x 

We can see from Eqs. (169) and (170) that the shape func- 
tions are built-up from an eight-node hexahedron representation 
(no midside nodes existing) . When we have no midside nodes we 
arrive at 
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In order to evaluate the terms in the 
the relationship 
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where J is given by Eq. (175). 


For an arbitrary element configuration the explicit integra- 
tion represented by Eq. (174) cannot be carried out. We employ 
a Gauss quadrature numerical integration scheme as an efficient 
procedure (Ref. 45). 


4.9.4 Geometric Stiffness Matrix 

There is no geometric stiffness matrix currently available 
with this element in PLANS. 
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4.9.5 


Initial Strain Stiffness Matrix 


The initial strain stiffness matrix is given by 

a a _i 


[k'] = 


‘’-r-r-1 


[B]^[C](det J)d^dCdTi 


(178) 


where J is given by Eq. (175) and [C] is defined in Sec- 
tion 4.9.8 for this element. To obtain Eq. (178) the initial 
strains are assumed to be constant throughout the element and 
evaluated at the element centroid. The plastic load vector is 
given by 

(Q) = [k*][€^} (179) 


4.9.6 Stresses and Strains 

The element stresses are found at the centroid from 

{a} = [CH(e^-€^} 

They are given in the principal material directions. The total 
strains are obtained from Eq. (176). The initial strains are 
due to both plastic and thermal effects and are the same as used 
in deriving the effective load. The elastic coefficient matrix, 
[C], is given in Section 4.9.8 for this element. 

4.9.7 Thermal Loads 

Thermal loads can be applied to an element by specifying 

the nodal temperature, T, and coefficients of expansion, a^. 

The thermal strains are e =e =e =a.T,€ =e = 

XX yy zz i xy yz 

e =0. These are treated as initial strains and the load 
xz 

vector is obtained from Eq. (179). 
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The nodal temperatures are input only at corner nodes and 
are averaged to obtain an element temperature for use in the 
thermal load vector. The may be different in each of the 

principal directions of orthotropy. 


4.9.8 Material Properties 

The elastic coefficient matrix, [C], corresponds to an 
orthotropic material and can be expressed as 


[C] = tQ]^[C^][Q] (180) 

where represents the reduced coefficient matrix along the 

principal axes of orthotropy and [Q] represents the transforma- 
tion into global coordinates. The terms of are as follows 
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The transformation is taken first as a rotation about the 
z-axis, then about the x-axis , then about the y-axis, i.e., 

Q = Q Q Q 

Three types of hardening laws may be specified for these 
elements. First, elas tic-ideally plastic behavior based on 
Hill's criterion for orthotropic materials may be input. This 
requires inputing the three yield stresses in the principal di- 
rections of orthotropy and the three yield stresses in shear 
with respect to these planes. The two other types of hardening 
may currently be used for initially isotropic material behavior. 
These are linear strain hardening and nonlinear strain hardening 
based on a Ramberg-Osgood representation of the stress-strain 
law. At the end of each half-cycle of loading the plastic mate- 
rial properties may be changed. 

4.9.9 Loadings 

Two different types of loading conditions can be applied. A 
consistent load can be applied across a face of an element. A 
concentrated load can be placed at a node. 

Surface Tractions — The consistenc load vector is obtained 
from the last term of Eq. (173) . For a consistent load we 
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represent the applied tractions in the same functional form as 
the generalized unknowns, i.e., 

(T) = [N][T^} (181) 


where (T^} represents the applied nodal tractions. Substitut- 
ing Eq. (181) into Eq. (173) and integrating in the local co- 
ordinate system we arrive at 


[F,] = 


[N]’^[N]Adq^dq^ 
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where t) and rj define the plane over which the traction is 
applied and 
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From Eq. (182) we see that the form of the load depends on 
whether midside nodes appear. For a face with no midside nodes 
the load varies linearly. When a midside node appears the load 
varies quadratically in the direction of the edge of the mid- 
side node. It is important to note that the input is the applied 
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tractions in Eq, (181), at nodes and not the total load. 

Concentrated Loads — A concentrated load can be placed at 
any node by simply specifying its magnitude (total load) and 
direction. Thus, in this case we are actually inputting F^. 

4.9.10 General Comments 

• Integration Order - A rectangular parallelepiped, in global 
coordinates, of eight nodes requires a 2x2x2 point Gauss 
integration while a rectangular parallelepiped of twenty nodes 
requires a 3x3x3 point Gauss integration. Care must, how- 
ever, be taken as the element diverges from a parallelepiped. 

An exact stiffness can no longer be calculated but a good ap- 
proximation can be found by Increasing the order of integration. 
Orders of integration which are too low lead to elements that 
are too flexible. Orders of integration that are too high cause 
roundoff problems. It is suggested that when a choice has to be 
made the order of integration be kept on the low side. This 
eliminates any roundoff errors and since element stiffnesses are 
generally too high the error in the stiffness matrix is in the 
right direction. In general, it is best to keep the element as 
rectangular as possible. 

• Orthotropic properties may only coincide with the global 
axis currently. 

• This element is available in the HEX module. 
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